MEDICAL IMAGING TECHNIQUES INCLUDING ADAPTIVE RECONSTRUCTION

Abstract

A computer-implemented method for adaptive image reconstruction can include: generating, by a processor, a preliminary reconstruction image using a preliminary reconstruction configuration; automatically adjusting the preliminary reconstruction configuration to an updated reconstruction configuration, by, at least: obtaining, by the processor, preliminary information from the preliminary reconstruction image; accessing, by the processor, a database of reconstruction configurations, the database providing a mapping of characteristics of images and objects in the images to reconstruction configurations; and performing, by the processor, a lookup operation to identify the updated reconstruction configuration based on the preliminary information; and generating, by the processor, a reconstruction image using the updated reconstruction configuration. Intermediate multifrequency images are generated during generating the preliminary reconstruction image and/or the reconstruction image and can be used to obtain preliminary information from the preliminary reconstruction image or reconstruction information from the reconstruction image, respectively.

Claims

1. A computer-implemented method for adaptive image reconstruction, comprising: generating, by a processor, a preliminary reconstruction image using a preliminary reconstruction configuration; automatically adjusting the preliminary reconstruction configuration to an updated reconstruction configuration, by, at least: obtaining, by the processor, preliminary information from the preliminary reconstruction image; accessing, by the processor, a database of reconstruction configurations, the database providing a mapping of characteristics of images and objects in the images to reconstruction configurations; and performing, by the processor, a lookup operation to identify the updated reconstruction configuration based on the preliminary information; and generating, by the processor, a reconstruction image using the updated reconstruction configuration.

2. The method of claim 1, wherein generating the reconstruction image using the updated reconstruction configuration comprises: applying an inverse scattering algorithm based on propagation through freespace; and applying a phase mask at each propagation step of the inverse scattering algorithm to give a total field comprising an incident field plus a scattered field.

3. The method of claim 2, wherein the inverse scattering algorithm is adapted to a particular coordinate system.

4. The method of claim 2, wherein a size of steps used in the inverse scattering algorithm for the generating of the reconstruction image is based on anisotropic pixel selection.

5. The method of claim 4, wherein at least two iterations of the inverse scattering algorithm use a different anisotropic pixel selection.

6. The method of claim 4, wherein the anisotropic pixel selection includes a distance in one direction between surfaces of the steps used in the inverse scattering algorithm.

7. The method of claim 6, wherein the generating of the reconstruction image is performed at multiple frequencies, wherein the distance in one direction between surfaces of the steps for the anisotropic pixel selection varies across the multiple frequencies.

8. The method of claim 7, wherein the distance in one direction between surfaces of the steps for the anisotropic pixel selection further varies across iterations for a particular frequency.

9. The method of claim 1, further comprising: receiving an indication of a coordinate system from a selection of a cartesian coordinate system and a curvilinear coordinate system; and filtering available reconstruction configurations at the database of reconstruction configurations according to the coordinate system before generating the preliminary reconstruction image using the preliminary reconstruction configuration.

10. The method of claim 1, wherein the adaptive imaging method is adapted to a curvilinear coordinate system, wherein available reconstruction configurations including for the preliminary reconstruction configuration and the updated reconstruction configuration are based on the curvilinear coordinate system.

11. The method of claim 10, further comprising: receiving data collected using a cylindrical array of transmitters and/or receivers; and transforming the data into cylindrical coordinates.

12. The method of claim 1, further comprising: receiving data collected using vertically limited diffraction-less beams.

13. The method of claim 12, wherein the preliminary reconstruction image and/or the reconstruction image is generated by applying a 2D reconstruction algorithm.

14. The method of claim 12, wherein the preliminary reconstruction image and/or the reconstruction image is generated by applying a 3D reconstruction algorithm.

15. The method of claim 12, wherein thickness of the diffraction-less beams are varied.

16. The method of claim 12, wherein the diffraction-less beams are rotationally limited.

17. The method of claim 1, further comprising receiving data captured by an imaging procedure using a multifrequency data collection method.

18. The method of claim 17, wherein the multifrequency data collection method includes transmitting a pulse containing multiple frequencies towards an object being imaged.

19. The method of claim 1, wherein the preliminary reconstruction configuration comprises a set of instructions, comprising forming a sequence of intermediate preliminary reconstruction images, wherein each intermediate preliminary reconstruction image has an increased resolution relative to its previous intermediate preliminary reconstruction image.

20. The method of claim 19, wherein the sequence of intermediate preliminary reconstruction images starts at a low frequency intermediate preliminary reconstruction image and proceeds in a stepwise fashion to higher frequency intermediate preliminary reconstruction images.

21. The method of claim 19, wherein the sequence of intermediate preliminary reconstruction images is carried out at frequencies in a non-monotonic progression.

22. The method of claim 1, wherein the updated reconstruction configuration comprises a set of instructions, comprising forming a sequence of intermediate reconstruction images, wherein each intermediate reconstruction image has an increased resolution relative to its previous intermediate reconstruction image.

23. The method of claim 22, wherein the sequence of intermediate reconstruction images starts at a low frequency intermediate reconstruction image and proceeds in a stepwise fashion to higher frequency intermediate reconstruction images.

24. The method of claim 22, wherein the sequence of intermediate reconstruction images is carried out at frequencies in a non-monotonic progression.

25. The method of claim 1, wherein the preliminary reconstruction configuration comprises a set of instructions, comprising: capturing an entire volume of a 3-dimensional image space; using a maximum distance between cross-sections of the 3-dimensional image space; using a limited range of frequencies to generate the preliminary reconstruction image; and using a limited number of iterations to generate the preliminary reconstruction image.

26. The method of claim 25, wherein the updated reconstruction configuration comprises a set of instructions, comprising: capturing a portion of the volume of the 3-dimensional image space; using a limited distance between cross-sections of the 3-dimensional image space; using an expanded range of frequencies to generate the reconstruction image; and using an increased number of iterations to generate the reconstruction image.

27. The method of claim 1, further comprising: detecting one or more errors in the preliminary reconstruction image in real time while generating the preliminary reconstruction image; stopping generating the preliminary reconstruction image; selecting a new preliminary reconstruction configuration from the database of reconstruction configurations based on the detected one or more errors in the preliminary reconstruction image; and restarting the method beginning with generating the preliminary reconstruction image using the new preliminary reconstruction configuration.

28. The method of claim 27, wherein generating the preliminary reconstruction image comprises assigning values to a plurality of voxels, wherein detecting one or more errors in the preliminary reconstruction image in real time while generating the preliminary reconstruction image comprises: evaluating values of voxels located near each other; and detecting an error if the values of voxels located near each other are above or below a voxel value threshold.

29. The method of claim 28, wherein the object being imaged is a patient's breast, wherein detecting an error if the values of voxels located near each other are above or below a voxel value threshold comprises detecting the presence of a silicone implant if the values of voxels located near each other are above a silicone implant voxel value threshold.

30. The method of claim 28, wherein the object being imaged is a patient's breast, wherein evaluating values of voxels located near each other comprises evaluating values of voxels near a patient's chest wall, wherein detecting an error if the values of voxels located near each other are above or below a voxel value threshold comprises detecting the presence of fatty tissue or non-biological material near the patient's chest wall if the values of voxels located near the patient's chest wall are above a chest wall voxel value threshold.

31. The method of claim 1, further comprising: detecting one or more errors in the reconstruction image in real time while generating the reconstruction image; and stopping generating the reconstruction image.

32. The method of claim 31, further comprising: selecting a new preliminary reconstruction configuration from the database of reconstruction configurations based on the detected one or more errors in the reconstruction image; and restarting the method beginning with generating the preliminary reconstruction image using the new preliminary reconstruction configuration.

33. The method of claim 31, further comprising: selecting a new reconstruction configuration from the database of reconstruction configurations based on the detected one or more errors in the reconstruction image; and restarting the method beginning with generating the reconstruction image using the new reconstruction configuration.

34. The method of claim 31, wherein generating the reconstruction image comprises large-scale optimization, wherein detecting one or more errors in the reconstruction image in real time while generating the reconstruction image comprises: evaluating values of a residual of the reconstruction image; and detecting a local minimum if the values of the residual are above a local minimum threshold value.

35. The method of claim 1, further comprising: obtaining, by the processor, reconstruction information from the reconstruction image; and outputting the reconstruction image and reconstruction information.

36. The method of claim 35, wherein the reconstruction image is a multifrequency image, wherein intermediate multifrequency images are generated during generating the reconstruction image, wherein obtaining, by the processor, reconstruction information from the reconstruction image comprises: calculating a frequency-independent wavenumber by removing a frequency dependent part of a wavenumber from a Helmholtz equation based model using information across frequencies of the intermediate multifrequency images.

37. The method of claim 36, wherein obtaining, by the processor, reconstruction information from the reconstruction image further comprises: estimating mass density using the frequency-independent wavenumber.

38. The method of claim 37, wherein obtaining, by the processor, reconstruction information from the reconstruction image further comprises: identifying objects using the estimated mass density.

39. The method of claim 35, wherein the reconstruction image is a multifrequency image, wherein intermediate multifrequency images are generated during generating the reconstruction image, wherein obtaining, by the processor, reconstruction information from the reconstruction image comprises: identifying frequency-dependent characteristics using values across frequencies of the intermediate multifrequency images, including attenuation and speed of sound.

40. The method of claim 39, wherein obtaining, by the processor, reconstruction information from the reconstruction image further comprises: estimating porosity based on the identified frequency-dependent characteristics.

41. The method of claim 1, wherein the preliminary reconstruction image is a multifrequency image, wherein intermediate multifrequency images are generated during generating the preliminary reconstruction image, wherein obtaining, by the processor, preliminary information from the preliminary reconstruction image comprises: calculating a frequency-independent wavenumber by removing a frequency dependent part of a wavenumber from a Helmholtz equation based model or paraxial approximation using information across frequencies of the intermediate multifrequency images.

42. The method of claim 35, further comprising: obtaining an attenuation image of the reconstruction image; generating, by the processor, a reflection image; performing a morphological operation with respect to the attenuation image to generate processed attenuation image; fusing the processed attenuation image with the reflection image to generate a final reflection image; and outputting the final reflection image.

43. The method of claim 42, wherein the morphological operation comprises erosion.

44. The method of claim 1, wherein the object being imaged is a patient's breast, wherein obtaining, by the processor, preliminary information from the preliminary reconstruction image comprises estimating mammographic density by: separating exterior voxels from breast voxels of the preliminary reconstruction image; segmenting high speed value breast voxels from other breast voxels of the preliminary reconstruction image to generate a first segmented image, wherein high speed value breast voxels are breast voxels having a speed value above a threshold; removing, from the first segmented image, high speed value breast voxels corresponding to skin tissue of the patient's breast to generate a second segmented image; and calculating mammographic density by determining a percentage of the high speed value breast voxels in the second segmented image.

45. The method of claim 44, wherein obtaining, by the processor, preliminary information from the preliminary reconstruction image further comprises obtaining a size of an object being imaged, wherein the database of reconstruction configurations provides a mapping of object size and mammographic density to reconstruction configurations.

46. The method of claim 1, wherein obtaining, by the processor, preliminary information from the preliminary reconstruction image comprises identifying a size of an object being imaged by: separating a volume of the object being imaged from an exterior volume; and quantifying the size of the object.

47. A system for adaptive image reconstruction, comprising: a processor; a storage system; a database of reconstruction configurations; a program stored at the storage system and comprising instructions that when executed by the processor direct the system to at least: generate a preliminary reconstruction image using a preliminary reconstruction configuration; automatically adjust the preliminary reconstruction configuration to an updated reconstruction configuration, by, at least: obtaining preliminary information from the preliminary reconstruction image; accessing the database of reconstruction configurations, the database providing a mapping of characteristics of images and objects in the images to reconstruction configurations; and performing a lookup operation to identify the updated reconstruction configuration based on the preliminary information; and generate a reconstruction image using the updated reconstruction configuration.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0009] FIG. 1 illustrates an ultrasound-based tomographic system that may be used to collect data used in the described image reconstruction and analysis techniques.

[0010] FIG. 2 illustrates a process flow diagram of a process that can be carried out by an acquisition control system.

[0011] FIG. 3 illustrates an example architecture for adaptive imaging.

[0012] FIGS. 4A and 4B illustrate an example computer-implemented method for adaptive image reconstruction.

[0013] FIG. 5 illustrates an example adaptive imaging process for generating transmission/attenuation image(s).

[0014] FIGS. 6A-6D illustrate results of a cylindrical coordinate system-based algorithm.

[0015] FIG. 7 shows reflection images of a breast, including the fusing of a processed attenuation image and reflection image.

[0016] FIGS. 8A and 8B illustrate light wave propagation for diffraction-less beams.

[0017] FIG. 9 shows a log plot of attenuation coefficient vs frequency from Bamber, 1986, Attenuation and absorption. Physical Principles of Medical Ultrasonics, ed. C R Hill.

[0018] FIGS. 10A-10D illustrate identification of tissue type through the power law variation of attenuation with frequency.

[0019] FIG. 11 illustrates propagation between a transmitter and receiver.

[0020] FIG. 12 shows a plot of Cs vs compressional speed.

[0021] FIGS. 13A-13D provide images showing application of porosity estimation using Biot as described herein.

[0022] FIG. 14 is schematic showing the Reflection coefficients relevant to the layered medium Green's function approach.

[0023] FIG. 15 illustrates a rectangular scattering matrix according to an embodiment of the present invention.

[0024] FIG. 16 illustrates a N-by-N scattering region according to an embodiment of the present invention.

[0025] FIG. 17 illustrates a N-by-N scattering coalesced region according to the present invention.

[0026] FIG. 18 shows a transducer coupling according to an embodiment of the present invention.

[0027] FIG. 19 shows the flowchart for the solution of the forward problem by means of the Parabolic FFT Marching method.

[0028] FIGS. 20A/B/C show the flowchart for the solution of the Jacobian of the forward problem for the Parabolic FFT Marching method.

[0029] FIGS. 21A/B/C show the flowchart for the action of the Hermitian conjugate of the Jacobian for the Parabolic FFT Marching method.

[0030] FIGS. 22A/B/C/D show the flowchart for the application of the conjugate gradient method applied to the generalization of the Propagation-Backpropagation method of Inverse Scattering.

[0031] FIG. 23 is the flowchart for the original Propagation-Backpropagation Method described in the paper A Propagation-Backpropagation Method for Ultrasound Tomography, Frank Natterer, which is included herein as reference.

[0032] FIGS. 24A/B/C/D/E show the flowchart for the brightness maximization based phase aberration correction algorithm.

[0033] FIG. 25 shows the basic Geometry for the paper A Propagation-Backpropagation Method for Ultrasound Tomography, Frank Natterer, which is included herein as reference. This is referred to as FIG. 1 in this paper. The square Q.sub.j is the square of sidelength 2 whose boundary is made up of .sub.j (side-scatter directions), .sub.j.sup. (Backscatter direction), and .sub.j.sup.+ (forward scatter direction). It encompasses the region 22, which contains the support of the object function . .sub.j is the direction of the incident wavefield propagation.

[0034] FIGS. 26A and 26B show an elongated rectangle Q.sub.j enclosing the region , with boundaries .sub.j (side-scatter directions), .sub.j.sup. (Backscatter direction), and .sub.j.sup.+ (forward scatter direction). .sub.j is the direction of the incident wavefield propagation (for FIG. 26A).

[0035] FIG. 27 shows the geometry of an acoustic transducer array illuminating an anatomical region through an aberrating layer of fat. A region of interest (ROI), selected by the user, is also shown. The transducer elements that contribute to the image in the ROI are denoted e.sub.m1 to e.sub.m2.

[0036] FIG. 28 shows the method for solving for the object function by means of a SQUARE nonlinear map, iteratively, in the scattered field domain.

[0037] FIG. 29 shows the iterative method for solving for the object function by means of a SQUARE nonlinear map in the object function () domain.

[0038] FIG. 30 shows the geometric set up for a typical calibration procedure.

[0039] FIG. 31 shows a comparison of the normalized total fields' (predicted and measured) magnitude versus receiver position.

[0040] FIG. 32 shows the comparison of the normalized total fields' phase. Here, the two fields are the predicted field using the starting values for the scattering parameters, and the measured field.

[0041] FIG. 33 shows the detailed model used for the receiver in the acoustic 2D case, in a typical calibration setup.

[0042] FIG. 34 shows the geometric setup used for the calibration of the receiver and transmitter.

[0043] FIG. 35 shows the scattering parameters used in the 2D acoustic calibration procedure typical of the calibration procedure used in inverse scattering, both for acoustic and EM modalities.

[0044] FIG. 36 illustrates a generic imaging algorithm, according to the present invention.

[0045] FIGS. 37A-37D detail an expanded imaging algorithm with respect to that shown in FIG. 36.

[0046] FIG. 38 shows the generic solution to the forward problem according to the present invention.

[0047] FIGS. 39A and 39B show the generic solution to the forward problem (or any square system) using a biconjugate gradient algorithm according to an embodiment of the present invention.

[0048] FIGS. 40A and 40B show the generic solution to the forward problem using the stabilized conjugate gradient algorithm, according to an embodiment of the present invention

[0049] FIG. 41 illustrates the application of the Lippmann-Schwinger Operator to the internal field estimate according to an embodiment of the present invention.

[0050] FIG. 42 illustrates the application of the Lippmann-Schwinger Operator in the presence of layering to the internal field estimate according to an embodiment of the present invention.

[0051] FIG. 43 illustrates the application of the Hermitian conjugate of the Lippmann-Schwinger operator.

[0052] FIG. 44 illustrates the scattering subroutine, according to an embodiment of the present invention.

[0053] FIG. 45 illustrates the propagation or transportation of fields from image space to detector position.

[0054] FIGS. 46A/B/C illustrate application of the Jacobian Matrix in the generic case (free spaceor no layering), according to an embodiment of the present invention

[0055] FIGS. 47A/B/C illustrates application of the Hermitian Jacobian Matrix in the generic case according to an embodiment of the present invention.

[0056] FIGS. 48A/B illustrates application of the Jacobian Matrix with correlations according to an embodiment of the present invention.

[0057] FIG. 49 illustrates the application of the second Lippmann-Schwinger operator in the generic case.

[0058] FIG. 50 illustrates the application of Hermitian conjugate of the second Lippmann Schwinger Operators

[0059] FIG. 51 the illustrates application the Hermitian conjugate of transportation of fields from image space to detector position.

[0060] FIGS. 52A/B illustrate application of the Hermitian conjugate of the Jacobian Matrix in the presence of layering according to an embodiment of the present invention.

[0061] FIG. 53 shows an example computing system through which X may be carried out.

DETAILED DESCRIPTION

[0062] Medical imaging techniques including adaptive reconstruction are provided. Adaptive reconstruction involves the adaptability of the algorithms and operations for imaging an object based on specified criteria and/or tests. Imaging techniques incorporating adaptive reconstruction can be referred to as adaptive imaging.

[0063] FIG. 1 illustrates an ultrasound-based tomographic system that may be used to collect data used in the described image reconstruction and analysis techniques. Referring to FIG. 1, imaging system 100 can include a transmitter 101, receiver 102, and transceiver 103. When operating in ultrasound frequencies, imaging system 100 can perform both reflection and transmission ultrasound methods to gather data. The reflection portion (e.g., transceiver 103) directs pulses of sound wave energy into tissues and receives the reflected energy from those pulseshence it is referred to as reflection ultrasound. Detection of the sound pulse energies on the opposite side of a tissue after it has passed through the tissue is referred to as transmission ultrasound.

[0064] The transmitter 101 and a receiver 102 are provided on opposite sides to enable the performing of transmission ultrasound. The transmitter 101 and the receiver 102 may be in the form of an array of transmitters and receivers. The transmitter array can emit broad-band plane pulses (e.g., 0.3-2 MHz) while the receiver array includes elements that digitize the time signal. The transceiver 103, which can include a set of reflection transducers, can be used to perform reflection measurements. The reflection transducers of the transceiver 103 can include transducers of varying focal lengths, providing a large depth of focus when combined. The focus can be fixed or variable in either or both directions. For example, in some cases, the focus may be in the vertical direction only, and the horizontal focus may be dynamic as controlled by a beamformer. In another embodiment, beamforming is included for the vertical direction.

[0065] The transmitter 101 and/or transceiver 103 can include a waveform generator that can produce repetitive narrow bandwidth signals or a wide bandwidth signal. The advantages of using a suitable wide bandwidth signal is that in one transmit-receive event, information at multiple frequencies of interest may be collected. Indeed, by using multiple frequencies, it will typically be possible to obtain more data for use in reconstructing the image. Accordingly, a computing system performing image reconstruction can receive data captured by an imaging procedure using a multifrequency data collection method, which includes transmitting a pulse containing multiple frequencies towards an object being imaged.

[0066] The lowest frequency used for the signal is selected such that the wavelength of the signal is not significantly smaller than the object to be imaged. In one embodiment, for example, the wavelength may be approximately twice the relevant characteristic length of the object (e.g., breast or other anatomical body part such as abdomen, extremities, etc.) or other smaller length. As a non-limiting example, a selected lowest frequency for imaging a breast may be 0.300 MHz. The highest frequency is selected such that the acoustic signal may be propagated through the object without being absorbed to such an extent as to render detection of the scattered signal impossible or impractical. Thus, depending upon the absorption properties and size of the object which is to be scanned, use of multiple frequencies within this range constraint will typically enhance the ability to more accurately reconstruct the image of the object. The use of multiple frequencies or signals containing many frequencies also has the advantage of obtaining data that may be used to accurately reconstruct frequency dependent material properties.

[0067] A coupling medium (e.g., in the form of a liquid or gel for matching of refractive indices) can be used between the object of interest and the imaging system 100. For example, a receptacle 110 can be provided to present a water (or other liquid or gel) bath in which a patient may rest at least the region of interest (e.g., the part 112 being imaged). Because the motion artifacts associated with patient movement can affect the image quality, mechanisms can be provided to facilitate retention and positioning of the breast or other body part. For example, in the case of breast tissue, an adhesive pad with a magnet can be placed near the nipple region of the breast and docked to a magnetized retention rod that gently holds the breast in a consistent position during the scan. As another example, a membrane over the bath between the breast and the liquid is used to hold the breast (and allow for alternative liquids in the bath).

[0068] 360 of data can be obtained through rotation of the system. The system (particularly arms containing the transmitter 101 and the receiver 102) may rotate 360 to acquire measurements from effectively all the angles (e.g., data sufficient to provide a 360 view even if not taken at every angle between 0 and) 360 and collect tomographic views of ultrasound wave data. The reflection transducer data can be collected with one or more horizontal reflection transducers of transceiver 103 that acquire data in steps or continuously as they rotate 360 along with the transmitter 101 and receiver 102. The transducers may also be tilted with a non-zero polar or azimuth angle.

[0069] In a specific implementation, the system rotates around the patient while both transmission and reflection information are captured. It is not necessary to acquire an entire 360 scan; images can be reconstructed with limited information. For example, a patient can lie prone with their breast pendent in a controlled temperature water bath (e.g., 31 C.) within the field of view of the transmitter 101, receiver 102, and transceiver 103 as the transmitter 101, receiver 102, and transceiver 103 rotate 360 around the patient. Then, in one example case 180 projections of ultrasound wave data may be obtained. In another example case, 200 to up to 360 projections or beyond of the ultrasound wave data may be obtained. After performing a rotation, an array chassis holding the transmitter 101, receiver 102, and transceiver 103 can be raised or lowered in a desired sequence to acquire data at another level. The sequence of levels for acquisition may be monotonic or non-monotonic, depending on implementation

[0070] Other detector configurations may be used. For example, additional detectors in a continuous or discontinuous ring or polygon configurations may be used. Of course, any configuration selected will have tradeoffs in speed and cost. In addition, in some cases, reflection arrays (the transducers for the reflection measurements) can do double-duty and perform independent transmission and receiver functions as well as reflection measurements. In another embodiment, the reflection array consists of a matrix of elements with multiple rows and columns.

[0071] An acquisition control system can be used to operate the various active components (e.g., the transducers) and can control their physical motion (when system 100 is arranged in a rotating configuration). An acquisition control system can automate a scan in response to a start signal from an operator. This automated acquisition process does not require operator interaction during the scanning procedure. Once the scan is complete, the acquisition control system (or other computing system having access to the data) can compute the reflection, speed of sound, and attenuation results from the collected data.

[0072] Accordingly, the data captured by system 100 can be used to reconstruct an image using inverse scattering techniques. Applications describing techniques for inverse scattering include U.S. Pat. Nos. 4,662,222; 5,339,282; 6,005,916; 5,588,032; 6,587,540; 6,636,584; 7,570,742; 7,684,846; 7,699,713; 7,771,360; 7,841,982; 8,246,543; and 8,366,617, which are hereby incorporated by reference in their entirety-except for that which is inconsistent with the apparatus and techniques described herein. In addition, the techniques described in Three-dimensional nonlinear inverse scattering: Quantitative transmission algorithms, refraction corrected reflection, scanner design and clinical results, Wiskin et al., Proceedings of Meetings on Acoustics, Vol. 19, 075001 (2013), are hereby incorporated by reference in their entirety.

[0073] Results can be provided to various computing systems including a viewing station and/or a picture archival and communication system (PACS). Thus, images can be automatically acquired, stored for processing, and available for physician review and interpretation at a review workstation.

[0074] FIG. 2 illustrates a process flow diagram of a process that can be carried out by an acquisition control system. Referring to FIG. 2, operation of imaging system 100 of FIG. 1 can be carried out using an acquisition control system performing process 200. Here, in response to receiving an indication to initiate automated scanning (e.g., from an operator or from some programmatic trigger), an acquisition control system can initialize and send a transmission wave from a specified angle about a patient (210), for example from one or more transmitters (such as transmitter 101). As the receiver(s) 102 sense the signal transmitting through the patient (220), raw transmission data 221 is captured. Then, spatially compounded extended depth of focus B mode scans, for example using transceiver(s) 103, are acquired (230) to obtain raw reflection data 231. Of course, in some cases, the B mode scans may be performed before the transmission ones. Additionally, in some embodiments, reflection transducers may have different focal lengths to extend the overall depth of focus within the imaging volume. In one embodiment the focal length is in the vertical direction. Of course, as mentioned above, the focus can be fixed or variable in either or both directions. In some cases, the transmitters (and transceivers) can be used to transmit a pulse containing multiple frequencies towards an object being imaged.

[0075] The acquisition control system determines whether the detectors are in the final position (240). For a rotating system, the acquisition control system can communicate with a motor control of the platform on which the active components are provided so that a current and/or next position of the platform is known and able to be actuated. For a fixed system, the acquisition control system determines the selection of the active arrays according to an activation program. In some implementations, the system may be a combination of rotating and fixed. For example, the system may be partially rotating and/or allow for vertical movement only. Accordingly, the detection of final position may be based on information provided by the motor control, position sensors, and/or a position program (e.g., using counter to determine whether appropriate number of scans have been carried out or following a predetermined pattern for activating transceivers). If the detectors are not in final position, the acquisition control system causes the array to be repositioned (250), for example, by causing the platform to rotate or by selecting an appropriate array of transceivers of a fixed platform configuration. After the array is repositioned, the transmission wave is sent (210) and received (220) so that the raw transmission data 221 is collected and the B mode scans can be acquired (230) for raw reflection data 231. This repeats until the detectors are determined to be in the final position.

[0076] FIG. 3 illustrates an example architecture for adaptive imaging. Referring to FIG. 3, an adaptive imaging architecture 300, which may be implemented on a computing system such as described with respect to FIG. 53, includes a reconstruction manager 310. Reconstruction manager 310 includes the instructions for controlling an adaptive imaging algorithm. For example, reconstruction manager 310 initiates appropriate preprocessing 320 of imaging data such as transmission data 312 and/or reflection data 314, directs generation of transmission/speed of sound/attenuation image(s) 330, directs generation of reflection image(s) 340, applies appropriate image evaluation processes 350 with respect to the generated images, and applies appropriate post processing 360 for the generated images. Reconstruction manager 310 enables automatic adjustments to the applied algorithms based on various aspects of the object being imaged in part based on the image evaluation processes 350 and the image itself such that each process can be adapted, allowing for dynamic selection of optimal parameters. As an illustrative example, as part of generating transmission/attenuation image(s) 330, the reconstruction manager provides an initial reconstruction configuration (e.g., from configuration directory 370) for generating an image and then uses information obtained from the image to update the reconstruction configuration so as to generate a better image. This process can include preliminary image reconstruction 332, preliminary reconstruction image evaluation 334, and image reconstruction. For example, reconstruction manager 310 can perform operations 410, 420, and 430 described with respect to FIGS. 4A and 4B. Within each reconstruction process (e.g., including preliminary image reconstruction 332 and image reconstruction 336), aspects of the reconstruction can be dynamic and adaptable based on various characteristics identified, for example, using specified tests. For example, reconstruction manager 310 can provide a particular reconstruction configuration to use for generating an image along with the particular tests to apply during the reconstruction that can cause changes to certain parameters and/or cause the process to stop and be restarted with the same or different reconstruction configuration. In the illustrated example, the preliminary image reconstruction 332 involves reconstruction configuration Config A 380 and test(s) 382; and the image reconstruction 33 involves reconstruction configuration Config B 390 and test(s) 392.

[0077] FIGS. 4A and 4B illustrate an example computer-implemented method for adaptive image reconstruction. Referring to FIG. 4A, a method 400 for adaptive image reconstruction can include generating (410), by a processor, a preliminary reconstruction image using a preliminary reconstruction configuration; automatically adjusting (420) the preliminary reconstruction configuration to an updated reconstruction configuration; and generating (430), by the processor, a reconstruction image using the updated reconstruction configuration.

[0078] Referring to FIG. 4B, automatically adjusting (420) the preliminary reconstruction configuration to an updated reconstruction configuration can include obtaining (422), by the processor, preliminary information from the preliminary reconstruction image; accessing (424), by the processor, a database of reconstruction configurations, the database providing a mapping of characteristics of images and objects in the images to reconstruction configurations; and performing (426), by the processor, a lookup operation to identify the updated reconstruction configuration based on the preliminary information.

[0079] Generating (430) a reconstruction image using the updated reconstruction configuration can include applying an inverse scattering algorithm based on propagation through freespace (which can be modeled as empty water); and applying a phase mask at each propagation step of the inverse scattering algorithm to give a total field comprising an incident field plus a scattered field. This method for generating the reconstruction image is applicable for generating (410) the preliminary reconstruction image as well. As will be described in more detail herein, the inverse scattering algorithm can be adapted to a particular coordinate system. In addition, various aspects of the inverse scattering algorithm can be adapted in accordance with the adaptive imaging processes described herein.

[0080] After the reconstruction image is generated, reconstruction information can be obtained from the reconstruction image; and the reconstruction image and reconstruction information can be output, for example, to a display.

[0081] Returning to FIG. 3, as an illustrative example implementing adaptive imaging architecture 300, the process may be started by command prompts from a technician. Certain preselected parameters may be selected to provide a starting point for operations and/or provide user input to the adaptive imaging processes. A preliminary configuration is selected from the configuration directory 370 by default, based on characteristics of the data, and/or based on input received from the technician. The configuration directory 370 contains all possible recipes. Recipes are sets of frequencies to be used in reconstruction of images from the transmission data. The recipe can include a range of frequencies and the intervals between each frequency. The recipe can also include the iterations at each frequency and each data level. Other aspects can be part of the recipes including parameters such as sections, anisotropic pixel selection, coordinate system, and stopping criteria.

[0082] Based at least in part on the recipes that may be selected for the preliminary configuration, preprocessing for the raw transmission data 312 is carried out. For example, preprocessing can include performing a Fourier transform. The Fourier transform is a mathematical operation that transforms the data from the time domain to the frequency domain. After the Fourier transform, magnitude is the y-axis and frequency is the x-axis for the data for operations carried out using the cartesian coordinate system. In another embodiment, the order of the data is changed to allow for quicker sequential access in memory. For a cylindrical or spherical coordinate system, additional preprocessing may be performed to transform coordinates to cylindrical or spherical. In some cases, the data is collected on a system that has the geometric shape corresponding to cylindrical or spherical or ellipsoidal coordinate surfaces.

[0083] After preprocessing, the generation of a preliminary image is carried out. The preliminary configuration enables an initial image to be generated quickly. The preliminary image is a relatively fast image, formed using relatively fewer frequencies, for relatively fewer iterations, and is generally blurry. In one embodiment, the image is formed using all levels of data. To form the preliminary image: all data levels may be used. Data levels refers to one particular level of data collection. In one embodiment, the array chassis (see e.g., discussion of system of FIG. 1) rotates through 360 degrees collecting transmission and/or reflection data at various azimuthal angles. In another embodiment, the rotational and vertical movement are combined so the concept of data level refers to one 360 degree rotation before the array changes azimuthal direction. In another embodiment, the arrays change direction before all 360 degrees are traversed. All data levels means every cross-section from the raw data, whether the cross-sections intercept the object or not. Taking all data levels will provide quantitative data regarding the size of the object, which can be included in the preliminary information obtained during the preliminary reconstruction image evaluation 334.

[0084] Some of the magnitude and frequency data from the preprocessed data are selected for each image, based on the recipe from preliminary configuration (e.g., Config A 380). In some cases, the preliminary reconstruction configuration includes a set of instructions involving capturing an entire volume of a 3-dimensional image space; using a maximum distance between cross-sections of the 3-dimensional image space; using a limited range of frequencies to generate the preliminary reconstruction image; and using a limited number of iterations to generate the preliminary reconstruction image. Other intermediate preliminary reconstruction images may be included.

[0085] As an illustrative example, the recipe can include first generating an initial attenuation image (time of flight), then generating a 2D image (e.g., based on the initial attenuation image), then generating a 3D coarse image (e.g., based on the 2D image), and then generating a 3D fine image (e.g., based on the 3D coarse image), where the 2D image is generated for each frequency of a set of frequencies selected between 0.3 MHz and 0.8 MHz, the 3D coarse image is generated for each frequency of a set of frequencies selected between 0.4 MHz and 0.8 MHZ, and the 3D fine image is generated for each frequency of a set of frequencies selected between 0.8 MHz and 1.3 MHz for example. Other embodiments can use different sets of frequencies, which may overlap for the three grids mentioned here. Note the coarse grid refers to voxels with generally larger dimensions and/or anisotropy.

[0086] Accordingly, multiple images are generated for the preliminary imageeach based on a different frequency signal/band. For each image, large-scale optimization can be performed (running inversions). In particular, as will be described in more detail herein, acoustic wave propagation through an object is simulated and the properties of the object are optimized to make the simulated propagation fit the observed propagation from the transmission data. The optimization problem can be solved using a stochastic gradient algorithm and the optimization results are used to assign values to pixels and form images. The pixels are squares on a 2D grid. Each 2D grid corresponds to a cross-section, or data level of the object being imaged. In some cases, as described in more detail herein, an adaptable anisotropic pixel may be used as a parameter in the reconstruction algorithm.

[0087] The reconstruction manager 310 provides the stopping criteria (e.g., included in test(s) 382) which in various embodiments may involve various data (e.g., residuals) associated with the minimization procedure, and subsets chosen for each iteration by the stochastic gradient algorithm and these characteristics can be dynamic for different iterations. The number of iterations for optimization may generally increase with each image formed (e.g., from the TOF, to the 2D image, to the coarse 3D image, to the fine 3D image), leading to higher quality images at the end of the sequence.

[0088] As part of the test(s) 382, detection of local minima can be performed early in the optimization process, for example, detection of local minima during the 3D coarse image generation where the residual is above a threshold can trigger the operations to cease and restart using a different recipe. This test of residual behavior may also be referred to as a process detecting whether a salt and pepper image is being formed-such an image is a pseudo-random juxtaposition of high and low values in place of physiologically correct estimates of the speed of sound. In such a case, the test may communicate this failure to the reconstruction manager 310, which can then obtain a different preliminary configuration, and restart the transmission image formation process.

[0089] Other tests 382 can also evaluate the object being imaged. For example, there may be a test to detect the presence of a silicone (or other) implant if there is a high residual or a high number of voxels which have a large number of pixels above a certain threshold (e.g., using pixel counts and speed of sound estimation). The test may determine the number of voxels having this value above the certain threshold is too large (above a certain number, for example). Tests can be performed to identify whether high speed voxels occur near the chest wall in the case of breast tissue (or other anomalous values that may be different than expected at a particular location) and this number is above a certain threshold. Indeed, tests 382 can include pixel (voxel) counts and speed of sound estimation, evaluation of gradients and residual magnitude, ratios of magnitudes of gradients and their behavior such as rate of decrease, and the like. Accordingly, tests can be included based on identified objects and based on the image formation process, which provide information used by the reconstruction manager 310 to select an updated reconstruction configuration (used to restart the preliminary image reconstruction 332 and/or generate the reconstruction image 336).

[0090] Indeed, in addition to stopping criteria and certain tests performed during the preliminary reconstruction, evaluation 334 of the resulting preliminary image is conducted. For example, for an image of a breast, quantitative measures such as percentage of fibroglandular tissue (sometimes referred to as mammographic density) or fibroglandular tissue volume or ratio (FGV, FGR)especially when skin volume is omitted) and the size of the breast can be obtained. These measures are then used to update the reconstruction configuration so as to best address the type of tissue and other object characteristics during the reconstruction algorithms. For example, the number of sections/data levels can be selected based on the size of the object/breast. Here, sections may refer to the number of data levels used in the reconstructionin one embodiment the top 24 mm (say 12 data levels at 2 mm separation) may be used for the reconstruction of the top section, the bottom 32 mm may be used for the reconstruction of the bottom part of the image and intermediate consecutive data levels may be used for the reconstruction of the middle part of the image. In addition, the frequency ranges may be adjusted based on the size, mammographic density, and results of test(s) 382.

[0091] In the following discussion, examples are not meant to be restrictive to a particular body part. Adaptive imaging is applicable to a variety of imaging applications including, but not limited to, the pediatric, orthopedic, and whole body imaging contexts. Tissue types vary with different applications (pediatric, etc.) and need not refer to breast imaging only. The presence of bone and air is a known complication that can be addressed using the described techniques. For example, the percentage of bone and air in a particular imaged volume can determine different recipes similar to that described in more detail below.

[0092] The updated reconstruction configuration (e.g., Config B 390) is then used to generate the reconstructed image 336. Compared to the preliminary image, the transmission image is a relatively slower image, using more frequencies, for more iterations, and is high-resolution. The process for forming the transmission image is similar to the process for the preliminary image. That is, acoustic wave propagation through an object is simulated and the properties of the object are optimized to make the simulated propagation fit the observed propagation from the transmission data. The optimization problem can be solved using a stochastic gradient algorithm and the optimization results are used to assign values to pixels and form images.

[0093] In some cases, the updated reconstruction configuration includes a set of instructions involving capturing a portion of the volume of the 3-dimensional image space; capturing a portion of the volume of the 3-dimensional image space; capturing a portion of the volume of the 3-dimensional image space; and using an increased number of iterations to generate the reconstruction image. As an illustrative example, the recipe for the updated reconstruction configuration can include generating a 3D coarse image, generating a first 3D fine image, and generating a second 3D fine image, where the 3D coarse image is generated for each frequency of a set of frequencies selected between 0.4 MHz and 0.8 MHz, the first 3D fine image is generated for each frequency of a set of frequencies selected between 0.8 MHz and 1.2 MHZ, and the second 3D fine image is generated for each frequency of a set of frequencies selected between 1.2 MHz and 1.3 MHz. In some cases, the 3D coarse image and/or 3D fine image generated during the preliminary image reconstruction 332 is used as a starting point for the 3D coarse image during the image reconstruction 336. As with the preliminary image, multiple images are generated based on a different frequency signal/band and the various adaptions and stopping criteria tests described with respect to the preliminary image reconstruction are applicable. More or fewer images and/or processes may be part of the updated reconstruction configuration (e.g., Config B 390) and tests 392. However, unlike for the preliminary image reconstruction, fewer than all data levels may be used. Here, the number of data levels/sections can be based on the size of the object.

[0094] FIG. 5 illustrates an example adaptive imaging process for generating transmission/attenuation image(s). Referring to FIG. 5, an example adaptive imaging process 500, which may be implemented as code executed by a computing system such as computing system 5300 of FIG. 53, can begin with an estimate of data properties 502. The estimate of data properties 502 may be received via a user input or be provided with the captured data (e.g., associated with a patient). From the estimate of data properties 502, an initial number of data levels (e.g., based on breast size) is determined for use by the fast reconstruction algorithm. An initial fast image is generated (504), for example using a coarse grid, 2D reconstruction, and/or time of flight algorithm. In some cases, the initial image is a transmission image. In some cases, a reflection image and/or data can be used to supplement the transmission image.

[0095] The generation of the initial fast image can be carried out such as described with respect to preliminary image reconstruction 332 of FIG. 3 and preliminary image reconstruction 410 of FIG. 4A. From the initial fast image, size of the breast can be estimated (506) and mammographic density estimated (508). The breast size estimation (506) can involve determining the size category for the breast, for example, using categories of small, medium, and large. The estimation of the mammographic density (508) involves determining the percentage of fibro-glandular tissue and determining the density category for the breast, for example, using categories of fatty, scattered/heterogenous, and dense. Based on the breast size and mammographic density, the recipe for inversion and number of sections can be determined (510). For example, the adaptive imaging process 500 can use a lookup table (LUT) or matrix, which in the case of breast can be such as follows.

Example LUT

TABLE-US-00001 size Density Small Medium Large Fatty Config 1 Config 2 Config 3 Scattered/heterogeneous Config 4 Config 5 Config 6 dense Config 7 Config 8 Config 9

[0096] In some cases, each configuration/recipe for the different combinations of breast tissue characteristics is different. In some cases, a same or similar recipe may be used for a subset of different combinations of breast tissue characteristics. Percentages of bone and air may also be part of the characteristics used in whole body/pediatric/partial body imaging. Similarly, for other anatomical regions of the body, a similar LUT could be created. For instance, for human body extremities and musculoskeletal tissue, it could be based on the expected ratio of fat vs muscle vs bone tissue.

[0097] Using the matrix/table one or more reconstruction configuration/recipes can be selected for use on each section of the breast. In the illustrated example, a coarse transmission image 512 is generated and a fine grid transmission image 514 is generated. These images can be generated such as described with respect to image reconstruction 336 of FIG. 3 and image reconstruction 430 of FIG. 4A.

Stopping Criteria and Tests

[0098] The reconstruction configurations can be selected and/or changed on the basis of various stopping criteria and tests (e.g., test(s) 382, 392). The criteria and tests can be based on image characteristics and can trigger the restart of an inversion process using a different set of parameters (e.g., frequencies, iterations, anisotropic pixels, or coordinate systems more suited to the data). The tests can include whether there is a silicone implant, whether there is a salt and pepper local minimum developing, and whether the image is going to be blurry.

[0099] Indeed, as mentioned above, as part of the test(s) 382, detection of local minima can be performed early in the optimization process, for example, detection of local minima during the 3D coarse image generation where the residual is above a threshold can trigger the operations to cease and restart using a different recipe. As another example in breast tissue, there is a test to detect the presence of a silicone implant if there is a high residual (e.g., using pixel counts and speed of sound estimation). Tests can be performed to identify whether abdominal fat is identified near the chest wall (or other anomalous values that may be different than expected at a particular location). Indeed, tests 382 can include pixel counts and speed of sound estimation, evaluation of gradients and residual magnitude, ratios of magnitudes of gradients and their behavior such as rate of decrease, and the like.

[0100] Thus, according to certain embodiments, an adaptive imaging method can thus include detecting one or more errors in the preliminary reconstruction image in real time while generating the preliminary reconstruction image; stopping generating the preliminary reconstruction image; selecting a new preliminary reconstruction configuration from the database of reconstruction configurations based on the detected one or more errors in the preliminary reconstruction image; and restarting the method beginning with generating the preliminary reconstruction image using the new preliminary reconstruction configuration.

[0101] In some cases, generating the preliminary reconstruction image comprises assigning values to a plurality of voxels, wherein detecting one or more errors in the preliminary reconstruction image in real time while generating the preliminary reconstruction image comprises: evaluating values of voxels located near each other; and detecting an error if the values of voxels located near each other are above or below a voxel value threshold. In an example implementation where the object being imaged is a patient's breast, detecting an error if the values of voxels located near each other are above or below a voxel value threshold can include detecting the presence of a silicone implant if the values of voxels located near each other are above a silicone implant voxel value threshold. Another test involving evaluating values of voxels located near each other includes evaluating values of voxels near a patient's chest wall, wherein detecting an error if the values of voxels located near each other are above or below a voxel value threshold comprises detecting the presence of fatty tissue or non-biological material near the patient's chest wall if the values of voxels located near the patient's chest wall are above a chest wall voxel value threshold.

[0102] Similarly, tests are applied during generating the updated reconstruction image. For example, the adaptive imaging method can include detecting one or more errors in the reconstruction image in real time while generating the reconstruction image; and stopping generating the reconstruction image. Similar to tests during the preliminary reconstruction, detecting one or more errors in the reconstruction image in real time while generating the reconstruction image can include evaluating values of a residual of the reconstruction image; and detecting a local minimum if the values of the residual are above a local minimum threshold value.

[0103] Once stopped, either a new preliminary reconstruction configuration can be selected (e.g., to start the process over at step 330) or a new updated reconstruction configuration can be selected (e.g., to restart process 350). For example, an adaptive imaging method can include selecting a new preliminary reconstruction configuration from the database of reconstruction configurations based on the detected one or more errors in the reconstruction image; and restarting the method beginning with generating the preliminary reconstruction image using the new preliminary reconstruction configuration. Of course, in other cases, the adaptive imaging method can include selecting a new reconstruction configuration from the database of reconstruction configurations based on the detected one or more errors in the reconstruction image; and restarting the method beginning with generating the reconstruction image using the new reconstruction configuration.

[0104] Although not described in detail, generation of reflection images can also include adaptive parameters and tests for the reflection image configurations can include, but are not limited to, whether there is a bright spot in the reflection image and whether there is a dark spot in the reflection image.

Adaptive Sectioning and Frequencies

[0105] The breast is a 3D object but has different characteristics relevant to wave propagation near the chest wall and in the sub-areolar regions. Based on the acquisition system, it is possible to include data redundancy in the vertical direction.

[0106] This means it is possible to collect data every 2 or 4 mm (which also could vary with the particular part of the volume of breast that is being imaged.) This vertical redundancy supports image quality in the presence of compromised SNR and data quality. By simultaneously inverting (imaging) from data at several different close levels, it is possible to increase the redundancy over not incorporating the information at the same time in inversion. It is noted that trying to invert from many levels simultaneously has an over-regularizing effect that creates a smoother image that may not be appropriate for the Reader or for diagnostic purposes. By separating the breast into sections, it is possible to optimize for the particular tissue/anatomy. In various embodiments, the breast can be divided into 2, 3, . . . . N sections.

Recipes

[0107] It is possible to determine a recipe for imaging the breast. A recipe includes the sequence of operations, including frequencies selected for use in the simulations. For example, a subset of frequencies within a range of frequencies corresponding to the frequencies used in the data acquisition can be selected (e.g., 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, . . . 1.35 MHz). The change between frequencies is determined empirically and is dynamic, it may change during the inversion process. In addition, the order in which frequencies are used is not required to be consecutive/monotonic That is, the frequency steps do not have to be monotonic increasing. In one embodiment the frequencies rise, then fall then continue to rise and this cycle can repeat. For example, a sequence can be 0.4, 0.425, 0.45, 0.475, 0.45, 0.425, 0.45, 0.475, 0.5, . . . .

[0108] In one embodiment a set of frequencies is chosen from 0.4, 0.425, . . . , 1.3 MHz, where the step/difference between frequencies is chosen here as 0.025 MHz. In another embodiment a subset of frequencies are imaged simultaneously. The spacing between the simultaneous frequencies being used can also be adaptive intra and inter-breast. That is, the spacing may change depending on the section. The code/program determines which set of frequencies and how many iterations to carry out at a given frequency to optimize image quality and also avoid local minima.

[0109] Accordingly, with reference to FIG. 3 and FIG. 4A, a preliminary reconstruction configuration can include a set of instructions, comprising forming a sequence of intermediate preliminary reconstruction images involving a sequence of intermediate preliminary reconstruction images, where each intermediate preliminary reconstruction image has an increased resolution relative to the previous intermediate preliminary reconstruction image. In some cases, the sequence of intermediate preliminary reconstruction images starts at a low frequency preliminary reconstruction image and proceeds in a stepwise fashion to higher frequency intermediate preliminary reconstruction images. In some cases, the sequence of intermediate preliminary reconstruction images is carried out at frequencies in a non-monotonic progression.

[0110] Similarly, an updated reconstruction configuration can include a set of instructions, comprising forming a sequence of intermediate reconstruction images involving a sequence of intermediate reconstruction images, where each intermediate reconstruction image has an increased resolution relative to the previous intermediate reconstruction image. In some cases, the sequence of intermediate reconstruction images starts at a low frequency reconstruction image and proceeds in a stepwise fashion to higher frequency intermediate reconstruction images. In some cases, the sequence of intermediate reconstruction images is carried out at frequencies in a non-monotonic progression. The sequence of frequencies for the updated reconstruction configuration may be the same or different than that used in the preliminary reconstruction configuration. In addition, for reconstruction configurations involving multiple levels of reconstruction (see e.g., examples provided for Config A 380 and Config B 390), where the intermediate images increase in resolution (e.g., from coarse to fine) the frequency sequence used in the reconstruction images may differ between the coarse reconstruction and the fine reconstruction.

Frequency Recipes

[0111] How the frequencies are chosen can impact whether local minima are encountered (see also stopping criteria). .sub.j, j=1, . . . , N.sub.k, where k=1, . . . , N.sub.freq. In one embodiment the frequencies are equally spaced. In another embodiment the frequencies are chosen so that delta lambda () is constant since the change in wavelength will determine whether cycle skipping takes place. In such an embodiment, to maintain a constant , the relationship between (frequency) and is not linear. That is, since

[00001]$f_j=-\frac{c_o}{{}_j^2}{}_j$

so that at a particular frequency j, the step size in frequency is

[00002]$f_j=-\frac{f_j}{{}_j}{}_j\mathrm{or}{f}_j=-\frac{f_j^2}{c_o}{}_j,$

that is the step size in frequency can in fact increase quadratically with frequency to keep the change in the wavelength the same. If the relative increase in wavelength is required to be constant than the formula

[00003]$f_j=-\left(\frac{{}_j}{{}_j}\right)f_j,$

indicates the step size increase linearly with frequency. The optimal step size may be determined empirically.

Adaptable Anisotropic Pixel

[0112] An anisotropic pixel selection can be carried out as part of an adaptive imaging process. That is, a size of steps used in an inverse scattering algorithm (e.g., propagation and backscattering) for generating a reconstruction image can be based on an anisotropic pixel selection, where the anisotropy is through the geometry (e.g., having different lengths in different directions). An anisotropic pixel selection includes a distance in one direction between surfaces of the steps used in the inverse scattering algorithm. The one direction may be in the direction of propagation, perpendicular to the direction of propagation or some other angle with respect to the direction of propagation. In some cases, at least two iterations of the inverse scattering algorithm use different anisotropic pixel selections. For example, the distance in one direction between surfaces of the steps used in the inverse scattering algorithm can change from one iteration to another iteration. As an illustrative example, a first distance may be used for five iterations and then a second distance is used for the next three iterations followed by a return to the first distance or a change to a third distance, which may be larger or smaller than the first and/or second distance.

[0113] The anisotropic pixel selection can vary between frequencies and iterations in any manner of combinations. In addition, to the distance in one direction that can vary, the direction in which that distance is adjusted may vary between frequencies and/or iterations.

Coordinate System

[0114] As mentioned above, the adaptive imaging process can be adapted to different coordinate systems. Available coordinate systems include, but are not limited to, rectangular, circular-cylinder, elliptic-cylinder, parabolic-cylinder, spherical, prolate spheroidal, oblate spheroidal, parabolic, conical, ellipsoidal, and paraboloidal. Indeed, available coordinate systems can be a cartesian coordinate system, a curvilinear coordinate system, or even a hybrid/expanding coordinate system. As long as an exact solution exists, it is possible to perform an inverse scattering algorithm based on propagation a short distance through free space (or empty water) via a closed or semi-closed form solution, for example, and application of a phase mask. The reconstruction configurations for the different coordinate systems utilize an appropriate phase mask. For example, cartesian coordinates have the plane wave, cylindrical coordinates have Bessel functions and exponentials, spherical coordinates have spherical Bessel (half integer order) and spherical harmonics (associated Legendre functions and exponentials), and expanding coordinates have narrow waist Gaussian functions multiplied by solutions to a paraxial approximation equation.

[0115] The ability to adapt the processes for a particular coordinate system such as curvilinear supports the receipt of data collected using a cylindrical array of transmitters and/or receivers. Such data can be transformed into cylindrical coordinates during preprocessing (e.g., operation 320 of FIG. 3) and subsequently analyzed in cylindrical coordinates.

[0116] In some cases, an indication of a particular coordinate system from a selection of a cartesian coordinate system and a curvilinear coordinate system can be received; and available reconstruction configurations at the database of reconstruction configurations can be filtered according to the particular coordinate system before generating the preliminary reconstruction image using the preliminary reconstruction configuration. In some cases, the system is adapted to a curvilinear coordinate system, wherein available reconstruction configurations including for the preliminary reconstruction configuration and the updated reconstruction configuration are based on the curvilinear coordinate system.

[0117] An example for spherical expanding coordinates via Gaussian beam is as follows.

[0118] Note the system involves propagation of energy primarily in the z direction. Here, the paraxial approximation equation can be obtained by factoring out the predominant wave propagating in the z direction:

[0119] Acoustic energy field: u(r)=u.sub.p(r)e.sup.iz.

[0120] Given the Helmholtz equation

[00004]$\frac{{}^{2}u}{x^2}+\frac{{}^{2}u}{y^2}+\frac{{}^{2}u}{z^2}+k^2\left(r\right)u=0$

[0121] We can make the assumption that the energy moves predominantly in the z direction and factor that out, as above, i.e.: u(r)=u.sub.p(r)e.sup.iz

[0122] This gives us

[00005]$\left(\frac{{}^2{u}_p}{x^2}+\frac{{}^2{u}_p}{y^2}-2i\frac{u_p}{z}-{}^2{u}_p+\frac{{}^2{u}_p}{z^2}+k^2{u}_p\right)e^{-iz}=0$

[0123] Finally, the variation in the acoustic field in the z direction is considered small compared to the variation that is factored out:

[00006]$.Math.\frac{u_p}{z}.Math..Math.\frac{{}^2{u}_p}{z^2}.Math.$

[0124] This gives the paraxial approximation:

[00007]$\frac{{}^2{u}_p}{x^2}+\frac{{}^2{u}_p}{y^2}-2i\frac{u_p}{z}+\left(k^2-{}^2\right)u_p=0$

[0125] Put =k for the following: (the energy is propagating at the appropriate frequency)

[0126] Gaussian Beam factor:

[0127] Note the Helmholtz equation does have the exact solution: (Green's function in 3D)

[00008]$u=\frac{e^{-\mathrm{ik}.Math.r-r^{}.Math.}}{4.Math.r-r^{}.Math.},$

[0128] Now using r=(0, 0, ib)

[0129] gives, with the Fresnel approximation:

[00009]$u=u_o\frac{{\mathrm{ibe}}^{-\mathrm{ik}\frac{x^2+y^2}{\left(z+\mathrm{ib}\right)}}}{\left(z+\mathrm{ib}\right)},$

[0130] Now put

[00010]$b=\frac{w_o^2}{2}.$

In this formulation this is a Gaussian beam with waist w.sub.o

[0131] So we thus have the diverging solution:

[00011]$u\left(x\right)=\frac{e^{-\mathrm{ik}\left(x^2+y^2\right)/q\left(z\right)}}{q\left(z\right)}$$\mathrm{Where}$$q\left(z\right)=z+i\frac{w_o^2}{}z+\mathrm{ib},b=\frac{2w_o^2}{2}=\frac{w_o^2}{2}$

[0132] Now factor out the Gaussian beam with w.sub.o waist size and use expanded coordinates for the remaining part of the wave, where for now is an arbitrary complex value:

[00012]$x^{}\frac{x}{z+\mathrm{ib}}$$y^{}\frac{y}{z+\mathrm{ib}}$$z_2^{}-z_1^{}=\frac{{}^2\left(z_2-z_1\right)}{\left(z_2+\mathrm{ib}\right)\left(z_1+\mathrm{ib}\right)}$$u\left(x\right)=\frac{v\left(x^{}\right)}{q\left(z\right)}{e}^{-\mathrm{ik}\left(x^2+y^2\right)/q\left(z\right)}$

[0133] The function variation (x) uses the expanding coordinates for reasons seen below:

[0134] We take a particular form of the diverging beam:

[0135] Now we assume that the Gaussian waist is infinitesimal, , so the wave is approximately spherical. Note the transformation:

[00013]$x^{}\frac{x}{z+i},y^{}\frac{y}{z+i},z_2^{}-z_1^{}=\frac{{}^2\left(z_2-z_1\right)}{\left(z_2+i\right)\left(z_1+i\right)}$

[0136] In the limit as e goes to zero and R ( real). will be constrained by the calculations below.

[00014]$x^{}=\underset{.fwdarw.0}{\lim }\frac{x}{z+i}=\frac{x}{z},y^{}\frac{y}{z},z_2^{}-z_1^{}=\frac{{}^2\left(z_2-z_1\right)}{z_2{z}_1}$

[0137] For calculation purposes we use:

[00015]$z^{}-z_1^{}=\frac{{}^2\left(z-z_1\right)}{\left(z+i\right)\left(z_1+i\right)}$

[0138] Use

[00016]$\frac{}{z}=\frac{x^{}}{z}\frac{}{x^{}}+\frac{y^{}}{z}\frac{}{y^{}}+\frac{z^{}}{z}\frac{}{z^{}}=\frac{x}{z^2}\frac{}{x^{}}+\frac{y}{z^2}\frac{}{y^{}}+\frac{{}^2}{z^2}\frac{}{z^{}}$

[0139] since

[00017]$\frac{z^{}}{z}=\frac{{}^2}{{\mathrm{zz}}_1}-\frac{{}^2}{z^2{z}_1}\left(z-z_1\right)=\frac{{}^2}{z^2}$$\mathrm{and}\frac{}{x}=\frac{x^{}}{x}\frac{}{x^{}}=\frac{}{z}\frac{}{x^{}}$$\frac{}{y}=\frac{y^{}}{y}\frac{}{y^{}}=\frac{}{z}\frac{}{y^{}}$

[0140] And so

[00018]$\frac{{}^2}{x^2}=\frac{}{z}\frac{x^{}}{x}\frac{{}^2}{x^2}={\left(\frac{}{z}\right)}^2\frac{{}^2}{x^2},$

with similar expression for y coordinate, and

[00019]$\frac{e^{-\mathrm{ik}\left(x^2+y^2\right)/q\left(z\right)}}{q\left(z\right)}$

solves the paraxial equation.

[0141] Now factor out the expanding wave and transform to the expanded coordinate system.

[00020]$u\left(x\right)=v\left(x^{}\right)w\left(x\right)=\frac{v\left(x^{}\right)}{q\left(z\right)}{e}^{-\mathrm{ik}\left(x^2+y^2\right)/q\left(z\right)},$

and we get.

[0142] Using,

[00021]$\frac{{}^{2}v\left(x^{}\right)w\left(x\right)}{x^2}=\frac{{}^{2}v\left(x^{}\right)}{x^2}w\left(x\right)+2\frac{v\left(x^{}\right)}{x}\frac{w\left(x\right)}{x}+\left(x^{}\right)\frac{{}^{2}w\left(x\right)}{x^2},$

and similar for y and chain rules above, so that:

[00022]$\frac{{}^{2}v}{x^2}+\frac{{}^{2}v}{y^2}-2\mathrm{ik}\frac{v}{z^{}}=0$

[0143] Propagation in the primed coordinates:

[0144] The total field is

[00023]$u\left(x\right)=\frac{v\left(x^{}\right)}{z}{e}^{-\mathrm{ik}\left(x^2+y^2\right)/z}$

[0145] The first fact satisfies the paraxial equation (by above)

[00024]$\frac{{}^{2}v}{x^2}+\frac{{}^{2}v}{y^2}-2\mathrm{ik}\frac{v}{z^{}}=0,$

so use transverse Fourier transform, F.sub.T:

[00025]$\stackrel{}{v}\left(k_{x^{}},k_{y^{}},z^{}\right)=\left(x^{},y^{},z^{}\right)e^{-k_T.Math.x^{}}{\mathrm{dx}}_T^{}=F_T\left(v\right)$

[0146] Equation is now

[00026]$\left(k_x^2+k_x^2\right)\hat{v}-2\mathrm{ik}\frac{\stackrel{}{v}}{z^{}}=0\mathrm{and}\mathrm{so}\frac{\stackrel{}{v}}{z^{}}=\frac{1}{2{\mathrm{ik}}_o}\left(k_x^2+k_x^2\right)\stackrel{}{v}$

has the exact solution:

[00027]$\stackrel{}{v}\left(k_{x^{}},k_{y^{}},z^{}\right)=e^{\frac{1}{2{\mathrm{ik}}_o}\left(k_x^2+k_x^2\right)\left(z^{}-z_o^{}\right)}A,$

constant A: that is:

[00028]$\stackrel{}{v}\left(k_{x^{}},k_{y^{}},z^{}\right)=e^{\frac{1}{2{\mathrm{ik}}_o}\left(k_x^2+k_x^2\right)\left(z^{}-z_o^{}\right)}\stackrel{}{v}\left(k_{x^{}},k_{y^{}},z_o^{}\right)$

[0147] Is the exact solution. Using z=(zz.sub.o)

[00029]$v\left(x^{},y^{},z^{}+z^{}\right)=F_T^{-1}\left(\stackrel{}{v}\left(k_{x^{}},k_{y^{}},z^{}+z^{}\right)\right),$

putting in all the operators:

[00030]$F_T^{-1}$

is the inverse transverse Fourier transform:

[00031]$\left(x^{},y^{},z^{}+z^{}\right)=F_T^{-1}\left(e^{\frac{1}{2{\mathrm{ik}}_o}\left(k_x^2+k_x^2\right)z^{}}\stackrel{}{v}\left(k_{x^{}},k_{y^{}},z^{}\right)\right)=F_T^{-1}\left(e^{\frac{1}{2{\mathrm{ik}}_o}\left(k_x^2+k_x^2\right)z^{}}{F}_T\left(v\left(x^{},y^{},z\right)\right)\right)$

[0148] In operator theoretic notation (read from right to left) this is just the transport in free space formula.

[00032]$\left(x^{},y^{},z^{}+z^{}\right)=F_T^{-1}\left[P_O\right]F_{T}v$

[0149] with free space propagator:

[00033]$P_o=e^{\frac{1}{2{\mathrm{ik}}_o}\left(k_x^2+k_x^2\right)z^{}}$

[0150] The total solution is:

[00034]$u\left(x\right)=\frac{v\left(x^{}\right)}{q\left(z\right)}{e}^{-\mathrm{ik}\left(x^2+y^2\right)/q\left(z\right)}$$q\left(z\right)=z+i\frac{w_o^2}{}=z+\mathrm{ib},b=\frac{2w_o^2}{2}=\frac{w_o^2}{2}$

[0151] Letting go to zero to get the used coordinate system:

[00035]$\underset{.fwdarw.0}{\lim }\frac{x}{z+i}=\frac{x}{z},y^{}=\frac{}{z},z_2^{}-z_1^{}=\frac{{}^2\left(z_2-z_1\right)}{z_2{z}_1}$

[0152] The coordinate transformation valid for z.sub.j1<z<z.sub.j

[0153] Gives for the original coordinates:

[00036]$y\frac{z}{}{y}^{},x\frac{z}{}{x}^{},;=z_i;R_i=z/z_i$$x^{}=\underset{.fwdarw.0}{\lim }\frac{x}{z+i}.fwdarw.x^{}=\frac{x}{R_i}$

[0154] z and x, y coordinates valid for z.sub.j1<z<z.sub.j

[0155] use =z.sub.j; R.sub.i=z/z.sub.j;

[0156] so now

[00037]$\begin{array}{c}{z}_j^{}-z_{j-1}^{}=\frac{{}^2\left(z_j-z_{j-1}\right)}{z_j{z}_{j-1}}=\frac{\left(z_j-z_{j-1}\right)}{R_{j-1}}\\ z_j^{}=\frac{z_j}{R_{j-1}};y^{}=\frac{y}{R_{i-1}};x^{}=\frac{x}{R_{i-1}}\\ z_j^{}=\frac{z_j}{R_{j-1}}=\frac{\left(z_j-z_{j-1}\right)}{R_{j-1}}\end{array}$

[0157] Propagate from z.sub.j1.fwdarw.z.sub.j in primed coordinates using the propagator derived above, then evaluate Gaussian part at

[00038]$y_{j+1}\frac{z}{}{y}_{j+1}^{}=R_i{y}_{j+1}^{};x=\frac{z}{}{x}^{};x=R_i{x}^{};$

and z=z.sub.j+z.sub.j=z.sub.j+1 (recall that z.sub.j1(z=z.sub.j1)=0 by construction so,

[00039]$\frac{e^{-ik\left(x_{j+1}^2+y_{j+1}^2\right)/2z_{j+1}}}{z_{j+1}},$

then use

[00040]$u\left(x\right)=v\left(x^{}\left(x_{j+1}\right)\right)\frac{e^{-ik\left(x_{j+1}^2+y_{j+1}^2\right)/2z_{j+1}}}{z_{j+1}}$

[0158] Then apply the phase mask for z.sub.j1.fwdarw.z.sub.j which is

[00041]$e^{{}_{z_{j-1}}^{z_j}{k}_o\mathrm{dz}},$

dz is the infinitesimal along the coordinate x, y path in x, y, z space:

[00042]$x=\frac{z}{}{x}^{};x=R_i{x}^{}\mathrm{and}y=\frac{z}{}{y}^{}=R_i{y}^{}$$u\left(x\right).fwdarw.u\left(x\right)e^{{}_{z_{j-1}}^{z_j}{k}_o\mathrm{dz}}$

[0159] In this way we march along for Z.sub.o, . . . , Z.sub.j1, Z.sub.j Notice the x, y coordinates do expand due to the R.sub.i1 and increasing as z increases.

[0160] FIGS. 6A-6D illustrate results of a cylindrical coordinate system-based algorithm.

[0161] Referring to FIG. 6A, a simulation model was created from a breast image and is shown along with the position of the cylindrical algorithm circular array and the circle defining the point source boundary condition.

[0162] The solution was first computed by solution of the integral equation. The simulation model was 1024 by 1024, wavelength pixels at 1.5 MHz (128 by 128 mm). The point source solution was computed with the Bi-Stab, FFT algorithm with a k-average preconditioner which converged in 40 steps. The total field on the receiver circle was then computed from the rectangular grid values using bi-quadratic interpolation.

[0163] The total field on the receiver circle was then computed via a cylindrical FFT-propagator, phase-mask algorithm.

[0164] FIG. 6B shows a comparison of a total field magnitude on the receiver circle; and FIG. 6C shows a magnified solution match. It can be seen that there is agreement of the cylindrical propagation algorithm with the simulation solution (referred to as IE solution) for all scattering angles.

[0165] FIG. 6D illustrates a cylindrical problem space with a radius of 100 mm and a height of 60 mm. Referring to FIG. 6D, the source is shown as a truncated, focused line source 30 mm high. The line source is also apodized with a Hamming window. The line source focal range=100 mm. The small, offset circle containing the source has a radius, a.sub.1=5 mm. The line source position in the small circle is x.sub.1=2.5 mm, y.sub.1=1.5 mm. This offset was done to ensure proper propagation of an asymmetric field.

[0166] The goal is now to evaluate the field on the 5 mm radius cylinder by Green's theorem and then to propagate the field onto the 100 mm radius cylinder using the cylindrical propagator algorithm. The Green's theorem calculation from the focused line source to the 5 mm cylinder is given by:

[00043]$f_1\left({}_1,z\right)={}_{-\frac{h_L}{2}}^{\frac{h_L}{2}}\frac{w\left(z^{}\right)e^{-{\mathrm{ik}}_0\sqrt{{\left(a_1\cos \left({}_1\right)-x_L\right)}^2+{\left(a_1\sin \left({}_1\right)-y_L\right)}^2+{\left(z-z^{}\right)}^2}-\sqrt{z^{{}^2}+f_L^2}}}{4\sqrt{{\left(a_1\cos \left({}_1\right)-x_L\right)}^2+{\left(a_1\sin \left({}_1\right)-y_L\right)}^2+{\left(z-z^{}\right)}^2}}{\mathrm{dz}}^{}$

[0167] where w(z) is a Hamming window, a.sub.1=5 mm, h.sub.L=30 mm, .sub.L=100 mm, x.sub.L=2.5 mm, y.sub.L=1.5 mm,

[00044]$k_0=\frac{2f}{c_0},f=1.5\mathrm{MHz},c_0=1.5\mathrm{mm}/\mathrm{microsecond}.$

[0168] Returning again to FIG. 3, after generating the transmission image, reconstruction manager 310 directs generation of reflection image(s) 340 from the reflection data 314. During the preprocessing 320, the intensity of acoustic reflection based on the time it takes the waves to travel back to the source can be obtained from the reflection data 314 (e.g., a Fourier transform may be applied such as described with respect to the transmission data 312). In addition, other preprocessing operations may be carried out such as coordinate system adjustments, etc. Any suitable reconstruction algorithm may be used to generate the reflection image. In some cases, the reflection image is corrected for refraction based on the transmission image. For example, more detail is added to interfaces within the object being imaged through use of the attenuation image (from the transmission data). The interface data can also be added after image formation. As an example implementation, an attenuation image can be obtained of the reconstruction image (e.g., generated in process 336), a reflection image can be generated (e.g., via suitable reconstruction algorithm), a morphological operation can be performed with respect to the attenuation image to generate a processed attenuation image, and the processed attenuation image can be fused with the reflection image to generate a final reflection image.

[0169] FIG. 7 shows reflection images of a breast, including the fusing of a processed attenuation image and reflection image. As shown in FIG. 7, the original reflection images correspond to the reflection image generated by any suitable reconstruction algorithm. Attenuation images are obtained and processed using a morphological operation. Examples of morphological operations that may be carried out on the attenuation images include, but are not limited to, erosion, dilation, opening, and closing. The processed attenuation images are then fused with the original reflection images to generate the final reflection images shown. As a comparison, the speed of sound (SOS) image is shown in the figure to illustrate that the information now visible in the final reflection image are not artifacts being added to the image. Advantageously, when the attenuation images highlight interface information, that interface information can improve the reflection image.

[0170] With the transmission image(s) and the reflection image(s) generated by processes 330 and 340, appropriate image evaluation processes 350 with respect to the generated images can be applied. The quantitative data obtained from the generated images can be similar to that obtained with respect to the preliminary image, but with more resolution/accuracy. For example, the quantitative data obtained from the generated images can be used for diagnostic and other purposes. Post processing 360 can be applied to the images to remove noise, perform decluttering, as well as apply desired metadata to make the images into a useful package for viewing and/or as training data for machine learning. As an example, the final reflection image(s) and transmission image(s) can be converted to DICOM format.

[0171] In some implementations, data for transmission data 312 and/or reflection data 314 can be acquired using a system such as described with respect to FIG. 1 and/or FIG. 2. Indeed, data captured by an imaging procedure using a multifrequency data collection method such as by transmitting a pulse containing multiple frequencies can be received.

[0172] In some cases, adaptive imaging architecture 300 includes or is in communication with a data resource of data collected using vertically limited diffraction-less beams. For such implementations, image reconstruction is carried out on the data collected using vertically limited diffraction-less beams. For example, in some of such implementations, the preliminary reconstruction image and/or the reconstruction image is generated by applying a 2D reconstruction algorithm. Through the adjusting of the reconstruction configuration as described herein, it is possible to reduce artifacts and speed up reconstruction.

[0173] In other of such implementations, the preliminary reconstruction image and/or the reconstruction image is generated by applying a 3D reconstruction algorithm. In some cases, the thickness of the diffraction-less beams are varied (and the 3D reconstruction can reflect such variance). In some cases, the diffraction-less beams are rotationally limited.

[0174] Creating diffraction-less beams to sonicate a breast (or other tissue) has the advantage that only second order scattering with be out of plane. The acoustic energy must scatter out of the plane, then scatter back into the plane of the receiver array.

[0175] FIGS. 8A and 8B illustrate light wave propagation for diffraction-less beams. FIG. 8A illustrates a cross-section of non-diverging diffraction-less beams and FIG. 8B illustrates a cross-section of a present acoustic total field. A system for creating diffraction-less beams can form beams similar to optical light sheets used in lightsheet microscopy. As illustrated in FIG. 8A, it can be seen that the beam remains narrow and can hold form for approximately 171, which would be 252 mm at 1 MHz (acoustic). These beams can be concatenated to create the horizontal sheets.

Formulation:

[0176] The well known expression generates a Gaussian profile Bessel beam

[00045]$\left(,z\right)-\frac{ikA}{2zQ\left(z\right)}{e}^{ik\left(z+\frac{{}^2}{2z}\right)}{J}_0\left(\frac{ik{k}_{}}{2zQ\left(z\right)}\right)e^{-\frac{1}{4Q\left(z\right)}\left(k_{}^2+\frac{k^2{}^2}{z^2}\right)}$

[0177] where Q(z)(qik/2z), k.sub. is transverse (x-y) wavenumber,

[0178] It is also possible to use superposition as these represent solutions to the wave equation. They are reasonable approximations when considering finite apertures.

[00046]${}_j\left(,z\right)-\frac{ik{A}_j}{2z{Q}_j\left(z\right)}{e}^{ik\left(z+\frac{{}^2}{2z}\right)}{J}_0\left(\frac{ik{k}_{}}{2z{Q}_j\left(z\right)}\right)e^{-\frac{1}{4Q\left(z\right)}\left(k_{}^2+\frac{k^2{}^2}{z^2}\right)}$

[0179] with Q.sub.j(z)(q.sub.jik/2z)

[0180] It is known these beams do not diffract or spread out as normal propagating wavefields do. This spreading requires the 3D algorithms to account for first order scattering from planes outside of the central plane. Only second order scattering effects occur with the confined beam. Using these Gaussian-Bessel beams (and superpositions of them) means that only second order scattering will occurscattering out of plane followed by scattering back into plane.

[0181] Such a superposition is:

[00047]$\left(,z\right){.Math.}_j{}_j\left(,z\right)$

[0182] This formulation can create a pseudo-plane wave with limited diffraction in the vertical direction (a sheet of acoustic energy). This is similar to light sheet microscopy fields.

[0183] It is also possible to use a time domain plane wave (pseudo plane wave).

[0184] Advantageously, the confinement of the acoustic signal means that out of plane scattering is second order and so a 2D algorithm of any kind will produce a much better image than with standard unconstrained ultrasound incident fields. These beams can be produced in the time domain as well, as constrained pulses and 2D plane wave pulses.

Evaluation Processes

[0185] As noted above with respect to operations for evaluating a preliminary reconstruction image 334 (including obtaining (422) preliminary information from the preliminary reconstruction image) and image evaluation 350, mammographic density can be measured (see also operation 508 of FIG. 5). In some cases, estimating mammographic density includes separating exterior voxels from breast voxels of the preliminary reconstruction image; segmenting high speed value breast voxels from other breast voxels of the preliminary reconstruction image to generate a first segmented image, wherein high speed value breast voxels are breast voxels having a speed value above a threshold; removing, from the first segmented image, high speed value breast voxels corresponding to skin tissue of the patient's breast to generate a second segmented image; and calculating mammographic density by determining a percentage of the high speed value breast voxels in the second segmented image.

[0186] Other information can also be obtained from the preliminary reconstruction image and/or final/updated reconstruction image. As also noted above, multifrequency images are collected. While generating a final multifrequency image, intermediate multifrequency images can be kept and used to obtain information from the preliminary reconstruction image and/or final/updated reconstruction image. As such, certain characteristics of the imaged object that vary with frequency can be used to extract useful information. While attenuation and speed of sound images are described in detail, the principle applies to images based on the backscatter (reflection) data, including with respect to attenuation coefficient slope estimate (ACS) and envelope statistics and derived parameters such as the effective scatterer diameter (ESD), effective acoustic concentration (EAC), and the parameter (effective scatterer density) from the HK distribution. Some of the estimates can be determined directly from the backscattered signal. As provided in detail herein certain parameters can be estimated based on tomographic reconstructed images using the transmitted signal.

Tissue Type Identification

[0187] It has been verified in the literature that certain mammalian tissue types behave in a particular manner with frequency. It is known that attenuation varies with frequency via the power law dependency. Advantageously, by using the Kramers-Kronig relationship that relates the frequency dependence of the attenuation to the speed of sound frequency dependence, it is possible to yield a corresponding power law for the speed of sound (and use the described reconstruction images to obtain useful information to identify tissue/lesions). From the power law relationship, it can be seen that the rate of change frequency for the SOS is less than for attenuation. In fact, the frequency dependence for speed can be shown to be

[00048]$\frac{1}{c\left(x,\right)}-\frac{1}{c\left(x,{}_o\right)}={}_o\left(x\right)\tan \left(\frac{a}{2}\right)\left({}^{a-1}-{}_o^{a-1}\right)$

when ()=.sub.o.sub.a (see below).

[0188] This is reflected in measured values as well, where the dependence of speed is hard to measure, whereas the dependence of attenuation is measured relatively easily. There is strong evidence from the literature that attenuation varies with frequency in a way that may help to characterize biological tissue.

[0189] In particular, the standard tissue model used for attenuation is

[00049]$\left(\right)={}_o{}^a,=2f$

[0190] which leads to the linear relationship:

[00050]$\log \left(x,\right)=\log {}_o\left(x\right)=a\log .$

[0191] This is acknowledged in the literature as being an approximation since the true model involves a summation of terms which represent relaxation processes at various frequencies. The phenomenological results of investigations indicate the power law to be generally satisfied.

[0192] FIG. 9 shows a log plot of attenuation coefficient vs frequency from Bamber, 1986, Attenuation and absorption. Physical Principles of Medical Ultrasonics, ed. C R Hill. The plot of FIG. 9 shows the power law behavior for most tissue types is a power law variance with frequency.

[0193] The tissue type (bone, skin, ducts, glands, fat, fibroglandular tissue, cancer, etc.) may affect the attenuation coefficient as well as the power law. Speed varies with frequency by virtue of the Kramers-Kronig relations. Note that these relations are non-local in their exact form. They involve integration over the entire real line in fact. There are, however, semi-local forms that can be used in an approximate manner. Note also that an argument based on operator series shows that for the special case of a power law variation of attenuation with frequency, i.e., ()=.sub.o.sub.a the speed of sound must also have variation with frequency given by the formula: (proved below)

[00051]$\frac{1}{c\left(x,\right)}-\frac{1}{c\left(x,{}_o\right)}={}_o\left(x\right)\tan \left(\frac{a}{2}\right)\left({}^{a-1}-{}_o^{a-1}\right).$

[0194] These formula show that the variation of speed with frequency will not necessarily be large since a1.3. This is borne out by experimental evidence as shown in FIGS. 10A-10D.

[0195] FIGS. 10A-10D illustrate identification of tissue type through the power law variation of attenuation with frequency. FIG. 10A shows a log-log plot of attenuation vs frequency (0.775, 0.8, 0.825, 0.85, 0.875, 0.9 MHz) for an image with a known cancer lesion. The region of inspection (ROI) is a 2 mm diameter centered in the known cancer lesion as determined by the speed of sound map. FIG. 10B shows a log-log plot of attenuation vs frequency for an image with fat. FIG. 10C shows panels of intermediate multifrequency images of a point in cancer and FIG. 10D shows panels of intermediate multifrequency images of a point in fat. Each image of the multifrequency images correspond to an image at a particular frequency (e.g., starting from the top with attenuation/SOS at 1.3 MHz and subsequent panels shown the attenuation images at sequentially lower frequencies to 0.400 MHz).

[0196] The calculation of these parameters over multiple frequencies can be carried out over a volume or region of interest to increase signal-to-noise (SNR) and stability. As can be seen, there is a different relationship between attenuation and frequency for a fat ROI as compared to a cancer ROI.

[0197] Using formula: log (x, )=log (x)= log

[0198] The coefficient here is a=10.2.

[0199] Power law result: Note these are empirical results subject to noise and are for illustrative purposes only. More accurate estimates can be obtained using a volume of interest in place of a single voxel.

TABLE-US-00002 Power law Fat 10.2 cancer 6.16

[0200] Advantageously, this strong separation through the power law parameters can be used to assign a cancer risk.

[0201] Accordingly, it is possible to obtain reconstruction information from a reconstruction image by identifying frequency-dependent characteristics using values across frequencies of the intermediate multifrequency images, including attenuation and speed of sound. For example, it is possible to calculate a frequency-independent wavenumber by removing a frequency dependent part of a wavenumber from a Helmholtz equation based model using information across frequencies of the multifrequency image. By quantitatively determining the spectral behavior of the speed of sound and attenuation images, it is possible to determine any correlation to the linear coefficients to the tissue type and/or abnormality. It is also possible to estimate mass density and porosity using such spectral behavior (e.g., the frequency-independent wavenumber). Indeed, as described herein, it is possible to estimate porosity based on identified frequency-dependent characteristics.

[0202] As an illustrative example of estimating mass density, once the attenuation is determined spatially (see e.g., examples in the imaging algorithms), mass density can be determined by first isolating the frequency independent wavenumber, solving a simple forward problem to estimate density distribution, and dividing out the density to yield the bulk modulus.

[0203] It should be noted that in the linear approximation (Born approximation) of the forward problem a simple calculation indicates that different parts of the data originate from either the monopole scattering from the bulk modulus or the dipole scattering from the density variations.

[0204] As illustrated, the data comes from different scattering potentials depending on where it is located: Note this is not an inversion scheme, rather it shows that there is some separation in the scattered data that may be exploited by a suitable scheme that has suitable regularization (see below for regularization schemes). Given a configuration such as shown in FIG. 11, which illustrates propagation between a transmitter and receiver. Note the Lippmann-Schwinger (LS) equation reads:

[00052]$f^i\left(r\right)=f\left(r\right)-k_o^2{}_{R^3}{}_{}\left(r^{}\right)f\left(r^{}\right)g_{k_o}\left(R\right){\mathrm{dr}}^{}+{}_{R^3}{}^{}.Math.\left({}_{}\left(r^{}\right)\left\{{}^{}f\left(r^{}\right)\right\}\right)g_{k_o}\left(R\right){\mathrm{dr}}^{}$

[0205] Directly from the differential equation. Integrate by parts and use Gauss theorem to remove total divergence:

[00053]$f^i\left(r\right)=f\left(r\right)-k_o^2{}_{}\left(r^{}\right)f\left(r^{}\right)g_{k_o}\left(R\right){\mathrm{dr}}^{}+{}_{}\left(r^{}\right)\left\{{}^{}f\left(r^{}\right)\right\}.Math.{}^{}{g}_{k_o}\left(R\right){\mathrm{dr}}^{}$

[0206] Using R=rr and =.sup.i(r)=e.sup.ik.sup.i.sup..Math.r (plane wave incident) and the standard approximation:

[00054]$g_{k_o}\left(R\right)==\frac{e^{-\mathrm{ikR}}}{4R}=\frac{e^{-\mathrm{ik}}.Math.r-r^{}.Math.}{4.Math.r-r^{}.Math.}\frac{e^{-ik\left(r-u_r.Math.r^{}\right)}}{4.Math.r.Math.}=\frac{e^{-ikr}}{4r}{e}^{-ik{u}_r.Math.r^{}}$

and

[0207] g.sub.k.sub.o (R)=iku.sub.re.sup.iku.sup.r.sup..Math.r for Green's function with approximation.

[0208] Define: k.sub.r=ku.sub.r vector pointing toward the receiver array element.

[0209] Substitute plane wave and Green's function approximation in the IE to get

[00055]$f^s\left(r\right)\frac{e^{-\mathrm{ikR}}}{4R}\left(k^2{}_{}\left(r^{}\right)e^{-i\left(k_r-k_i\right).Math.r^{}}{\mathrm{dr}}^{}+k_i.Math.k_r{}_{}\left(r^{}\right)\left\{e^{-i\left(k_r-k_i\right).Math.r^{}}\right\}{\mathrm{dr}}^{}\right)$

[0210] Inner product is cos of angle between them times k.sup.2 so:

[00056]$f^s\left(r\right)\frac{e^{-\mathrm{ikr}}}{4r}{k}^2\left({}_{}\left(r^{}\right)e^{-i\left(k_r-k_i\right).Math.r^{}}{\mathrm{dr}}^{}+\cos {}_{}\left(r^{}\right)\left\{e^{-i\left(k_r-k_i\right).Math.r^{}}\right\}{\mathrm{dr}}^{}\right)$

[0211] These are Fourier transforms and give expression for scattered field.

[00057]$f^s\left(r\right)\frac{e^{-\mathrm{ikr}}}{4r}{k}^2\left({}_{}\left(k_r-k_i\right)+\cos {}_{}\left(k_r-k_i\right)\right)$

where is the angle between k.sub.r, k.sub.i.

[0212] This gives the relation in the linear approximation between scattered data and the object functions for density and bulk modulus. To complete the separation, the following correlations between mass density and speed can be utilized, particularly for the speed values of interest for mammalian tissue.

[00058]$K=\left(1706.507c^2-1611.737\right)\mathrm{MPa}$$K=893c^3-349c^2=c^2\left(893c-349\right)\mathrm{MPa}$

[0213] The above correlations are particularly useful as a regularizing term in density reconstruction and can be used during the preliminary image reconstruction. For example, it is possible to estimate the bulk modulus and shear modulus for breast and mammalian tissue from these correlations. Note:

[0214] K=c.sup.2c.sup.3 which indicates that the cube of the speed image can be considered a stiffness image.

[0215] Bulk modulus and shear wave estimation is possible from the speed of sound images (e.g., as reconstruction information from the reconstruction images). For example, the shear modulus turns out to be more sensitive in some sense to the presence of a hard lesion, since a cancer can have a shear speed 30 times greater than normal tissue.

[0216] Accepting (based on

[00059]$c_s=\sqrt{\frac{}{}})$$\frac{c}{c_s}=\sqrt{\frac{K}{}}$

[0217] gives (see below)

[00060]$\frac{K}{9000}$

[0218] so that

[00061]$\frac{c^3}{9000}$

[0219] Accordingly, it can be seen that a shear moduli estimation is possible.

[0220] In particular, the estimate of based on the shear wave speed of tissue being approximately 5 m/sec is

[00062]$=\frac{K}{{\left(c/c_s\right)}^2}\frac{K}{{\left(300\right)}^2}=\frac{K}{9*10^4}$

[0221] The Bulk modulus formula:

[00063]$K=893c^3=c^2\left(893c-349\right),c\mathrm{in}\mathrm{mm}/\sec ,K\mathrm{in}\mathrm{MPa}.$

[0222] Using this formula, it is possible to show stiffness directly in the images.

[0223] Here stiffness is bulk modulus. Thus it is possible to show a quantitative representation of stiffness that others are showing qualitatively. Note that it is possible to determine the Young's modulus as 3, the shear modulus presumably using the perturbation formula from (E is Young's modulus):

[00064]$E=\frac{3+2}{+}3\left(1+2/3\right)\left(1-\frac{}{}+.Math.\right)$

[0224] The perturbation expansion to first order in being E3.

[0225] Note also, as illustrated by the plot of FIG. 12, the shear wave speed greatly magnifies the differences in speed between compressional wave speed 1560 to 1610. This is the range for differentiating between fibroadenomas and cancer.

[0226] Using the formula:

[00065]$=\frac{K}{{\left(c/c_s\right)}^2},\mathrm{where}K=893c^3-349c^2=c^2\left(893c-349\right)\mathrm{MPa}$

Mass Density Estimation

[0227] Accordingly, it is possible to utilize the multifrequency images in attenuation and remove the frequency dependent part of the equivalent wavenumber that results from the square root transformation utilized to justify neglect of density. Once the frequency independent of the equivalent wavenumber is determined, there remains a nonlinear second order differential equation that can be transformed in two stages into a linear second order equation (essentially a Helmholtz equation) with known wavenumber and boundary conditions known (it is an elliptic equation).

[0228] The problem can be solved once (since it is not an inverse problem), and the spatial distribution of the density can be retrieved.

[0229] Once the density is known it is possible to estimate the bulk modulus by virtue of the formula:

[00066]$c=\sqrt{\frac{K}{}}$

[0230] The density dependent wave equation is known to be

[00067]${}^{2}p+\frac{}{c_o^2}p=.Math.\left({}_{}p\right)-\frac{}{c_o^2}{}_{}p$

where .sub.K(KK.sub.o)/K.sub.o, (.sub.o)/ are the gamma object functions for compressibility and density respectively. Pressure is p, frequency is =2f and background speed of sound is c.sub.o. Wavenumber

[00068]$k_o^2{}^2/c_o^2.$

Using

[00069]$p^{}=\frac{p}{\sqrt{}}$

the equation becomes .sup.2p(x)+k.sup.2(x)p(x)=0 with pseudo-wavenumber given by:

[00070]$k^2\left(x\right)=k^2\left(x\right)-\frac{3}{4{}^2}{\left(\right)}^2+\frac{1}{2}{}^2.$

Note the true wavenumber

[00071]$k\left(x\right)\frac{}{c\left(x\right)}+i\left(x\right).$

As noted the attenuation generally follows a power law with exponent of about 1.2 or so for mammalian tissue. Accordingly, it is possible to remove the frequency dependent part. Kramer's Kronig relations can be exactly integrated in this case and the power law for speed of sound 0.2 k.sub.ind (x) is the frequency independent part that remains when the frequency dependent part is subtracted out.

[0231] Because multiple frequency images are generated, there is access to this information. We use the identity

[00072]$\frac{1}{}{}^2=.Math.\left(\frac{1}{}\right)+\frac{1}{{}^2}.Math.:$

(proved by direct calculation).

[0232] The equation now reads:

[00073]$k^2\left(x\right)-k^2\left(x\right)k_{\mathrm{ind}}^2=-\frac{1}{4}\ln .Math.\ln +\frac{1}{2}.$

( ln ), To solve this equation we use the transformation: To solve this equation, we first transform to a new y variable:

[00074]$\frac{1}{2}{}_i\ln =-{}_i\ln y.$

The equation now reads

[00075]$-\frac{{}^{2}y}{y}=k_{\mathrm{ind}}^2\left(x\right)$

which can be rewritten as

[00076]${}^{2}y\left(x\right)+k_{\mathrm{ind}}^2\left(x\right)y\left(x\right)=0.$

[0233] The method can include generating images at the multiple frequencies; estimating the attenuation at multiple frequencies, removing the frequency dependent part of the wavenumber, leaving k.sub.ind(x), solving

[00077]${}^{2}y\left(x\right)+k_{\mathrm{ind}}^2\left(x\right)y\left(x\right)=0,$

the forward Helmholtz problem for y (suitable BCs), and estimating density as

[00078]$=\frac{e^c}{y^2}=\frac{K}{y^2},$

for suitably determined constant, K determined by the known density of water at the boundary and in the water bath.
Biotporous media

[0234] It is noted that the integral equation formulation for the Biot Theory in the Acoustic approximation is provided in the imaging algorithm section (see the vector Lippmann-Schwinger equationFredholm of Second Kind) as well as the paraxial approximation to the Helmholtz.

[0235] Expanding upon the imaging algorithm provided below, a conversion to the paraxial form is provided. An explicit method for determination of the porosity parameters based on the spectral determination of speed of sound and attenuation is included.

[0236] The full elastic Biot model is shown as follows:

[00079]$\left(\begin{array}{cc}{\hat{}}_{11}& {\hat{}}_{12}\\ {\hat{}}_{21}& {\hat{}}_{22}\end{array}\right)\left(\begin{array}{c}{}_t^{2}u\\ {}_t^{2}U\end{array}\right)=\left(\begin{array}{cc}P& Q\\ Q& R\end{array}\right)\left(\begin{array}{c}\left(.Math.u\right)\\ \left(.Math.U\right)\end{array}\right)-\left(\begin{array}{c}N\left(u\right)\\ 0\end{array}\right)$

[0237] This is reduced via homogenization theory to effective parameters (see below for definition of the parameters):

[00080]$U_s+\left(+\right)\mathrm{grad}\mathrm{div}\left(U_s\right)=-{}^2\left({}_{11}{U}_s+{}_{12}{U}_l\right)+i\left(U_s-U_l\right)$$n/\mathrm{grad}\mathrm{div}\left[\left(/n-1\right)U_s-U_l\right]=-{}^2\left({}_{21}{U}_s+{}_{22}{U}_l\right)-i\left(U_s-U_l\right),$$\begin{array}{ccc}\mathrm{strain}& \mathrm{inertial}& \mathrm{viscous}\\ \mathrm{energy}& \mathrm{terms}& \mathrm{dissipation}\end{array}$$\mathrm{where}$$H=H_1+{\mathrm{iH}}_2=1/K\left(\right)$

with

[00081]${}_{11}=\left(1-n\right){}_s+\left({\mathrm{nH}}_2/-{}_l\right)$$\begin{array}{cc}{}_{12}=n\left({}_l-H_2/\right)& ={\mathrm{nH}}_1\end{array}$${}_{21}=\left({}_l-{\mathrm{nH}}_2/\right)$${}_{22}={\mathrm{nH}}_2/.$

[0238] This is a full elastic representation and the details are in Boutin et al. (C. Boutin, G. Bonnet, and P. Y. Bard, Green functions and associated sources in infinite and stratified poroelastic media, Geophysical Journal International, vol. 90, pp. 521-550, 1987).

[0239] FIGS. 13A-13D provide images showing application of porosity estimation using Biot as described herein.

[0240] As shown in FIG. 13A, segmentation of marrow in tibia based on the SOS and constrained by ellipsoid gives average value 1402.2 m/s which agrees with literature values. FIG. 13B shows bone segmentation (trabecular) in tibia. FIG. 13C shows segmentation of trabecular bone in coronal and sagittal views in femur.

[0241] The Table 1 below gives the quantitative accuracy of the bone marrow and speed of sound of the Biot slow wave.

TABLE-US-00003 TABLE 1 comparison of segmented vs literature values for marrow and Biot slow wave human bone marrow Biot slow wave SOS Literature values QTUS measured values ~1410 m/s 1402, 1389 ~1470 m/s [2, 3] 1472, 1466 m/s

Homogenization in Biot Context:

[0242] Note that often in mathematical analysis of wave or diffusion phenomena in disordered or periodically structured media, homogenization is used. When the characteristic size of the periodicity for example is size <<1 compared with the O(1) size of the macrostructure, analytic expressions relating coefficients at the macroscopic scale based on the microstructure and the concomitant coefficients can be formulated. Furthermore, as the ratio .fwdarw.0, characteristic homogenization parameters are created. These macroscopic parameters are useful for studying wave and diffusion parameters in microstructured media. In particular in porous media:

[0243] Using the notation of Boutin et al., the following model for wave propagation in porous media results:

[00082]$\left(+2\right){}^2{P}_s+{}^2{P}_s-{}_0{}^{2}P=-.Math.F\left(x\right){}^{2}P-P-{}_0{P}_s=-V\left(x\right)$

[0244] Where P.sub.sU.sub.s is the divergence of particle motion in the solid matrix, and .sub.0, , are parameter's defined in the following way:

[00083]$=\left(1-n\right){}_s+{}_l\left[n+{}_l{}^{2}K\left(\right)/\mathrm{.Math.}\right]=K\left(\right)/\mathrm{.Math.}{}_o=+{}_l{}^{2}K\left(\right)/\mathrm{.Math.}=+{}_l{}^2$

[0245] Where, =1K.sub.b/K.sub.s, with K.sub.b=3+2/3, the bulk modulus, and l, m are the Lame' parameters in standard linear elastic theory. Also =n/K.sub.s+n/K.sub., and: [0246] is the frequency of the interrogating wave [0247] .sub.l is the density of the liquid phase [0248] .sub.s is the density of the solid phase [0249] n is a saturation parameter [0250] K() the generalized, explicitly frequency-dependent Darcy coefficient introduced via homogenization theory

[0251] In the above, an approximate value based on homogenization theory for straight ducts can be used:

[00084]$K\left(\right)=\left(\mathrm{ik}\right)/\left(v{}^{*}\right)J_2\left(-8i{}^{*}\right)/J_0\left(-8i{}^{*}\right);k=n{a}^2/8;{}^{*}=k/\mathrm{nv}$

with J.sub.0, J.sub.2 being the Bessel functions, a the radius of the ducts, and the kinematic viscosity of the fluid

[0252] This is an acoustic approximation to the full elastic equations but presents an opportunity for an exact solution to the Green's function equation:

BG.sub.Biot+I.sub.22(x)=0, where I.sub.22 is the 2D identity and

G.sub.Biot: R.sup.2,3.fwdarw.R.sup.2

[0253] Is a matrix Biot Green's operator. B is the operator matrix

[00085]$B=\left[\begin{array}{cc}{}^2-& -{}_0\\ -{}_0{}^2& {}^2+{}^2\end{array}\right]$

[0254] Using the adjoint of this matrix (i.e., the signed matrix of cofactors, not the linear adjoint), the following solution is obtained for the Biot Green's operator.

[00086]$G_{\mathrm{Biot}}=\frac{1}{4}\frac{{}_0}{\left({}_2^2-{}_1^2\right)}\left\{\left(w.Math.d\right)\left[\begin{array}{cc}& 0\\ {}_0& \end{array}\right]+\left(w.Math.v\right)\left[\begin{array}{cc}\stackrel{\_}{}{}^2& {}_0\\ 0& -\stackrel{\_}{}\end{array}\right]\right\}$

[0255] This is used in the Generalized acoustic Biot Lippmann-Schwinger equation:

[00087]$r_0\left(x\right)=r\left(x\right)-\left[\begin{array}{cc}-{\stackrel{\_}{b}}_1& 0\\ 0& {\stackrel{\_}{b}}_2\end{array}\right]G_{2x2}\left(x-\right)D^{}\left(\right)r\left(\right)d^3{D}^{}\left(x\right)=\left[\begin{array}{cc}{}_1& 0\\ 0& {}_2\end{array}\right]$

[0256] Contains the object functions that contain the tissue characteristics. Also where r: R.sup.2,3.fwdarw.R.sup.2, and s: R.sup.2,3.fwdarw.R.sup.2 are vector-valued functions on R.sup.3 or R.sup.2 given by

[00088]$r\left(x\right)=\left[\begin{array}{c}P\left(x\right)\\ P_s\left(x\right)\end{array}\right],s=\left[\begin{array}{c}V\left(x\right)\\ .Math.F\left(x\right)\end{array}\right]$

where w: R.sup.3.fwdarw.R.sup.2 is the 2D vector defined by:

[00089]$w\begin{array}{c}\frac{e^{-i{}_1\left|x\right|}}{.Math.x.Math.}\\ \frac{e^{-i{}_2\left|x\right|}}{.Math.x.Math.}\end{array}$

[0257] Where it is seen that the d.sub.i are effective wavenumbers in porous media and the vectors

[00090]$d=\left(\begin{array}{c}-{}_1^2\\ {}_2^2\end{array}\right),v=\left(\begin{array}{c}1\\ -1\end{array}\right)$

Are so defined.

[0258] The .sub.i are the effective wavenumbers.

[0259] Using this Green's function the forward problem is well defined.

[0260] Furthermore, the Jacobian calculation and adjoint of the Jacobian calculation proceeds as described herein (see e.g., step 3850 of FIG. 38 and SCIENTIFIC BACKGROUND ON SQUARING AN OVERDETERMINED SYSTEM TO APPLY BiSTAB)

[0261] These operations respectively provide a step length and the gradient direction.

[0262] This allows for the full implementation of the inversion algorithm.

[0263] In one embodiment the imaged parameters are: Darcy coefficient and Porosity. FIG. 13D shows an attenuation image of human (ex vivo) cadaver knee. The patella, tibia, fibula and femur are present, the tibiofemoral space is clear. This anomalous attenuation is related to the Biot generalized DArcy coefficient.

[0264] Other embodiments can include more or fewer parameters in addition to the standard wave speed and attenuation. The other parameters that are not imaged via this disclosure can be determined from literature values.

[0265] These calculations can be carried out on CUDA or similar AMD cards/processors for speed.

[0266] Acoustic theory appropriate for multiple organs

[0267] Note that although orthopedic images are shown here, the porosity and density parameters apply to lungs, kidneys, liver, and other organs or abdominal imaging.

[0268] Paraxial approximation for Biot wave theory:

[0269] Recall that pressure in the solid is: P.sub.sU.sub.s P.sub.sU.sub.s The inhomogeneous Biot equation from above is:

[00091]$\left(\begin{array}{cc}\left[\left(+2\right){}^2+{}^2\right]& -{}_o{}^2\\ -{}_o& {}^2-\end{array}\right)\left(\begin{array}{c}{P}_s\\ P\end{array}\right)=\left(\begin{array}{c}-F\\ -V\end{array}\right)$

[0270] Rewriting this as:

[00092]$\left(\begin{array}{cc}\left(+2\right)& -{}_o\\ 0& \end{array}\right){}^2+\left(\begin{array}{cc}{}^2& 0\\ -{}_o& -\end{array}\right)\left(\begin{array}{c}{P}_s\\ P\end{array}\right)=\left(\begin{array}{c}-F\\ -V\end{array}\right)$

[0271] And using the matrix definitions:

[00093]$A\left(\begin{array}{cc}\left(+2\right)& -{}_o\\ 0& \end{array}\right),B\left(\begin{array}{cc}{}^2& 0\\ -{}_o& -\end{array}\right)$

[0272] Gives the following form for the Biot acoustic approximation.

[00094]$\left(A{}^2+B\right)\left(\begin{array}{c}{P}_s\\ P\end{array}\right)=0$

[0273] This can be solved using the Biot Acoustic Green's function developed above. This leads to a generalized Lippmann-Schwinger equation that is rigorous and leads to an inversion algorithm by itself and the concomitant step length and gradient direction calculations.

[0274] As can be seen from the images, it is possible to detect the slow Biot nonstandard compression wave (P wave) and have segmented the trabecular bone region interior to the bone and obtained a value commensurate with the predicted speed of sound for this slow Biot wave, yielding clinically valuable effective parameters.

Paraxial Approximation:

[0275] This vector Helmholtz equation can be approximately factored into

[00095]$\left(\begin{array}{cc}\left(+2\right)& -{}_o\\ 0& \end{array}\right){}^2+\left(\begin{array}{cc}{}^2& 0\\ -{}_o& -\end{array}\right)\left(\sqrt{A}+i\sqrt{B}\right)\left(\sqrt{A}-i\sqrt{B}\right)$

[0276] Where

[0277] =(n)/K.sub.s+n/K.sub. is a volume averaged averaged Bulk modulus

[0278] And a is a relative difference modulus, and a.sub.o is a frequency dependent relative difference modulus with the effective Darcy coefficient K:

[00096]$=1-K_b/K_s,{}_o+{}_l{}^{2}K\left(\right)/i$

[0279] The matrix square root can be found in the standard way for A and B:

[0280] This leads to the parabolic PDE (paraxial approximation);

[00097]$\left(\sqrt{A}-i\sqrt{B}\right)P0$

[0281] By direct calculation then:

[00098]$\sqrt{B}=\left(\begin{array}{cc}\sqrt{}& 0\\ \frac{{}_o\left(i\sqrt{}-\sqrt{}\right)}{{}^2+}& i\sqrt{}\end{array}\right),\mathrm{and}$$\sqrt{A}=\left(\begin{array}{cc}\sqrt{+2}& \frac{{}_o\left(\sqrt{}-\sqrt{+2}\right)}{\left(+2\right)-}\\ 0& \sqrt{}\end{array}\right)$