WO2024095273A1

WO2024095273A1 - High quality high resolution 3-d medical imaging using ultrasound

Info

Publication number: WO2024095273A1
Application number: PCT/IL2023/051140
Authority: WO
Inventors: Anat Canning
Original assignee: Ultravista Ltd
Current assignee: Ultravista Ltd
Priority date: 2022-11-06
Filing date: 2023-11-06
Publication date: 2024-05-10
Anticipated expiration: 2025-05-06

Abstract

A method of obtaining a medical ultrasound image of a field of view of a subject's body, comprising: a) placing ultrasound source and receiver transducers on the body, surrounding the field of view; b) recording ultrasound signals transmitted by a source transducer through the body and received by a receiver transducer; c) using a velocity model to find travel times between grid points and source and receiver transducers; d) constructing gathers as a function of grid point and a parameter related to ray path length, using pre-stack depth Kirchhoff migration; e) correcting the velocity model using information from the gathers about how ray path length affects travel time; f) constructing a pre- stack attribute of the field of view as a function of grid point, from the gathers; and g) reconstructing an image of the field of view from the pre- stack attribute as a function of grid point.

Description

HIGH QUALITY HIGH RESOLUTION 3-D MEDICAL IMAGING USING ULTRASOUND

RELATED APPLICATIONS

This application claims the benefit of priority under 35 USC § 119(e) of U.S. Provisional Patent Application No. 63/423,007, filed on November 6 2022, the contents of which are incorporated herein by reference.

FIELD AND BACKGROUND OF THE INVENTION

The present invention, in some embodiments thereof, relates to an ultrasound imaging system and method, and, more particularly, but not exclusively, to a medical ultrasound system and method for generating an ultrasound image of a field of view of the interior of a human or animal body.

Sound waves are used for seismic imaging of geological regions of interest, for example for oil exploration and other mineralogical exploration. Various methods have been developed for constructing seismic images of underground regions from received acoustic echo data, including various migration methods, such as Kirchhoff migration, as well as full waveform inversion (FWI) methods. Reviews of migration methods for seismic imaging are given, for example, by Biondo L. Biondi, 3D Seismic Imaging, Investigations in Geophysics Series No. 14, Society of Exploration Geophysicists, Tulsa, Oklahoma, 2006; by A. J. Berkhout, Seismic Migration: Imaging of Acoustic Energy by Wavefield Extrapolation, A. Theoretical Aspects , second edition, Elsevier, 1982; and by Etienne Robein, Velocities, Time-imaging and Depth-imaging: Principles and Methods, third edition, European Association of Geoscientists and Engineers, 2010. Reviews of FWI methods in seismic imaging are given, for example, by J. Virieux and S. Operto, “An overview of fullwaveform inversion in exploration geophysics,” Geophysics 74, WCC127-WCC152 (2009), and by Virieux et al, “An introduction to full waveform inversion,” Encyclopedia of Exploration Geophysics, “Seismic inversion,” pages Rl-1 to Rl-30, SEG, 2017, DOI: 10.1190/1.9781560803027.entry6.

Taner et al, “Complex seismic trace analysis,” Geophysics 44, 1041-1063 (1979), explains the application of complex trace analysis to seismic data, and its usefulness in geologic interpretation. The conventional seismic trace can be viewed as the real component of a complex trace which can be uniquely calculated under usual conditions. The complex trace permits the unique separation of envelope amplitude and phase information.

Herbert W. Swan, “Velocities from amplitude variations with offset,” Geophysics 66, 1735- 1743 (2001), describes how the sensitivity of amplitude variation with offset (AVO) to normal moveout (NMO) velocity errors can be exploited to obtain more accurate velocities from a fully automated procedure.

Anat Canning and Alex Malkin, “Automatic anisotropic velocity analysis for full azimuth gathers using AV AZ,” SEG Houston 2009 International Exposition and Annual Meeting, p. 201- 205, describes an automatic velocity picking method for deriving high resolution velocities in seismic imaging. The method can be used to pick isotropic and anisotropic velocity fields simultaneously on full azimuth 3D angle gathers. The resulting velocity field is very detailed both laterally and vertically.

Ultrasound medical imaging has generally been done using a hand-held array of ultrasound transducers, for example a linear array, that sends a focused pulse of ultrasound waves into soft body tissue, and receives echoes reflected and scattered from the tissue. An image of a 2-D field of view is generally displayed in real time, using delay-and-sum beam forming, which uses the echo time to determine the depth of a feature, assuming a uniform speed of sound in the tissue. The plane of the 2-D field of view can be changed by manually manipulating the array. Reviews of such conventional ultrasound imaging methods are given by Jprgcn Arendt Jensen, “Medical Ultrasound Imaging,” Progress in Biophysics and Molecular Biology 93, 153-165 (2007), Elsevier; and by Peter N. T. Wells and Hai-Dong Liang, “Medical ultrasound: imaging of soft tissue strain and elasticity,” J. Royal Soc. Interface 8, 1521-1549 (2011).

Theis and Bonomi, “Seismic imaging method for medical ultrasound systems,” Physical Review Applied 14, 034020 (2020), describe a proposed ultrasound-imaging strategy, an application of concepts developed in the context of seismic prospecting, in contrast to conventional medical ultrasound systems that use delay-and-sum beamforming with an implicit assumption of constant velocity media. They use a 144-element linear transducer array, with 0.245mm interelement spacing. Echoes are sensed in a synthetic aperture configuration, with the transducer elements sequentially emitting a nearly spherical wave front that covers the whole region of interest.

Renaud et al, “Real-time Kirchhoff migration for ultrasound imaging of the bone cortex,” SEG International Exposition and 88^th Annual Meeting, p.4797-4801, 2018

(DOI: 10.1190/segam2018-2996779.1), explore the feasibility of real-time ultrasound imaging of the cortex of a long bone. They show that it is possible in vivo with a conventional clinical ultrasound probe and a fully programmable ultrasound scanner. The ultrasound images are reconstructed using Kirchhoff migration and taking into account the elastic anisotropy of cortical bone tissue. An image rate of 1 image per second was achieved. Published U.S. Patent Application 2019/0328355 to Oscar Calderon Agudo, and the related published PCT application WO 2018/108820, describe a non-invasive method of generating image data of intra-cranial tissue using ultrasound energy that is transmitted across a head of a subject through the skull of the subject. An ultrasound data set is derived from a measurement of ultrasound waveforms generated by a source of ultrasound energy, transmitted through the skull and intracranial cavity, and/or reflected by the inner or outer boundaries of the skull, detected by a plurality of receivers located at an opposing side of a region within the intra-cranial cavity. The ultrasound sources and receivers may be included in a helmet or other wearable device. A predicted data set is obtained by modeling data from a starting model of the acoustic properties of the area investigated, such as wave speed velocity. To build the model, an MRI image was used in combination with tables of known acoustic velocities of various brain and skull tissues. The waveforms in the observed data set are then analysed by comparing them to the waveforms in the predicted data set in order to recover a model of at least one acoustic property of the body. This can be done using the full-waveform inversion (FWI) method by minimizing the least squares norm between observed and predicted data.

Guasch, Calderon Agudo, et al, “Full-waveform inversion imaging of the human brain,” Nature Partner Journals: Digital Medicine (2020) 28, describe computer simulations to demonstrate that full-waveform inversion, a computational technique originally developed in geophysics, is able to generate accurate three-dimensional images of the brain with sub-millimeter resolution. The authors describe a prototype helmet containing 1024 transducers held rigidly in a 3-D printed network, customized to provide an accurate fit to an individual subject, as well as a portable clinical device in which the transducers move radially and contact the patient via a sonographic gel.

Published U.S. Patent Application 2021/0080573 to Bachmann et al, describes a system and method for producing ultrasound images based on Full Waveform Inversion (FWI). The system produces 3D maps of physical parameters that control wave propagation, such as shear and compressional wave speeds, mass density, attenuation, Poisson ratio, bulk and shear moduli, and impedance, for medical imaging and non-destructive testing. The system includes a transmitting element, which may have a precisely known and controlled location within the system, and a receiving element with a precisely known and controlled location within the system.

Robins et al, “Deep-Leaming-Driven Full-Waveform Inversion for Ultrasound Breast Imaging,” Sensors 2021, 21, 4570, describes using ultrasound computed tomography (USCT) for breast imaging, with image reconstruction using full-waveform inversion (FWI), an iterative data- fitting method that was originally developed for seismic exploration. Data for a breast phantom immersed in water was acquired using dual 96-element P4-1 cardiac probes, rotating the receiver probe relative to a transmitting element to simulate a fixed ring array of ultrasound transducers. A convolutional neural network (CNN) was used to extrapolate the received signals to frequencies below the range of sensitivity of the transducers, in order to allow the FWI to reconstruct an image. Some other papers on USCT are cited in <en(dot)wikipedia(dot)org/wiki/Ultrasound_computer_tomography>, downloaded 21 September 2022.

Muller et al, “Reverse Time Migration: A Seismic Imaging Technique Applied to Ultrasonic Data,” International Journal of Geophysics, July 2012, doi: 10.1155/2012/128465, describes fast and efficient multi-channel ultrasonic echo measurement systems introduced for practical applications in non-destructive testing of concrete elements (NDT-CE). The acquired multi-offset data are processed and imaged by the synthetic aperture focusing technique (SAFT), which is closely related to the seismic Kirchhoff migration method. Some limitations of SAFT are discussed.

An article in the MIT News, downloaded on October 12, 2022 from <news(dot)mit(dot)edu(slash)ultrasound-stickers-0728>, describes a rigid square ultrasound array designed by MIT engineers, 2 cm across, that can be stuck to the skin using a stretchy adhesive layer surrounding a layer of an elastic and stretchy solid hydrogel that transmits ultrasound. The array remains in place for up to 48 hours, and produces images of the inside of the body, including major blood vessels, the stomach, and the heart. The article describes earlier designs for patches with stretchable flexible ultrasound arrays, but such arrays are said to produce low-resolution images, due in part to the change in relative position of the transducers when the array stretches, and they cannot image deep organs.

International patent application publication WO2022/123464 to Axelrod and Chaviv describes a wearable sensing garment for sensing ultrasound waves from a body part. The sensing garment includes a mesh of fabric adapted to be worn on at least a body part of a wearer, and a plurality of ultrasound sensors disposed on the mesh fabric, wherein each of the plurality of ultrasound sensors are adapted to detect ultrasound waves returned from at least a body part of the wearer and further adapted to transmit at least a portion of the detected ultrasound waves to a controller. The ultrasound sensors are adapted to collect signals emitted by a plurality of ultrasound emitters randomly or deliberately embedded into the wearable garment. The controller is adapted to determine the positions of the sensors and mapping them with respect to the body parts based on returned soundwaves collected by the sensors, using adaptations of full wave inversion (FWI). The ultrasound sensing garment is designed to provide an imaging solution that transmits and receives ultrasound signals that are then processed to generate a high resolution, three-dimensional (3D), image of a scanned body part.

Published patent application US 2023/0157673 to Axelrod and Dafni describe a system and method for generating an ultrasound image of a moving medium. An initial model of elastic properties and a morphing function is generated to account for the changes of the medium motion over time, and continuously refined during the inversion process until its accuracy exceeds a predetermined threshold. The motions within the medium are small motions such as bowel motion, expansion of blood vessels, raising of the rib cage and the like. The disclosed embodiments apply to a variety of methods that internally use full wave propagation, e.g. FWI and its variants, Least Squares Migration, one- or two-way wave equation migration, Reverse Time Migration, other seismic migration variants, and more. Figure 1 depicts an example system for performing full waveform inversion (FWI) in a moving medium. The system comprises an ultrasound sensing garment, and further details are said to be described by Axelrod and Chaviv, cited above.

US Patent 11,559,285 to Tai et al describes a medical ultrasound (US) imaging system that includes a US probe and a processor. The US probe includes an array of transducers arranged in a reflection geometry, the probe configured to emit US waves and to receive reflected ultrasound waves that are reflected from a body portion of a patient. The processor is configured to generate an image of the body portion of the patient by applying an inverse model to the emitted and reflected US waves. In one embodiment, the image comprises at least one physiological tissue parameter generated by applying full waveform inversion (FWI) to the emitted and reflected soundwaves. In another embodiment, the processor is configured to an image of the body portion by applying reverse time migration (RTM) model to the emitted and reflected ultrasound waves.

Published patent application US 2015/0374334 to Klock et al describe methods and systems for imaging dense anatomical structures. It has generally been understood that ultrasound energy cannot penetrate bone and other more dense bodily structures. As a result, prior art reflection ultrasound imaging has been limited to soft tissues that are readily penetrable. However, the inventors have discovered that ultrasound energy traversing through bone can be detected using inverse scattering techniques. The inventors have discovered that transmission wave calculations, when altered to include shear wave moduli, can effectively model the wave energy conversion from compressional wave to shear wave and back again, thereby converting the weak transmitted signal into a usable signal for inverse scattering. An ultrasound scanning system for imaging a body part of a subject, which can be a joint and/or a bone, can include an imaging chamber operable to contain a liquid transmission medium and to receive a bone or joint of the subject into the transmission medium, and at least one transducer array disposed in the chamber. Published patent application US 2020/0253580 to Lenox et al describe Quantitative Transmission Ultrasound (QTUS) techniques and systems for the detection and determination of lesions in dense tissues. In certain implementations, QTUS uses ultrasound energy to image and characterize breast tissue. A QTUS system performs both transmission and reflection ultrasound methods to gather data. A receptacle is 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 being imaged). A transmitter and a receiver are provided on opposite sides to perform transmission ultrasound. A series of reflection transducers are also provided to perform reflection measurements.

US Patent 6,385,474 to Rather et al relates to a method and apparatus for the construction and/or use of multidimensional fields that can be used for high resolution detection and characterization of features within objects. The multidimensional field is collected by an array of radiation detectors that substantially surround the object. The detected radiation is produced by an array of radiation sources and is subsequently scattered, reflected, transmitted by the object under study and any features within the object under study. In particular embodiments of the invention, the radiation that is used is ultrasonic radiation and the object under study is human or animal tissue or an organ. In this case, the invention permits the detection and identification of cancer. In one embodiment, the arrays comprise large numbers of sources and detectors, with over 5000 detectors. With a sufficient number of such sources and detectors, the invention provides for construction of a three-dimensional rendering of numerous physical quantities to describe the object. If the object under study is a human breast, a modified hemispherical arrangement of sources will provide large coverage. In the case of study of the human breast, it is preferable to have additional coverage beyond a hemisphere, such as the geodesic dome configuration shown in Fig. 1(b). The detectors are also placed so as to cover a large solid angle, preferably greater than 2TI steradians.

EP 1578274 to Duric et al describes a sensor system comprising ultrasonic transmission and receiving elements, and a mechanism for transferring both transmitted and received acoustic energy to the tissue being studied. One embodiment uses a pair of acoustic paddles used with breast imaging, both upper and lower paddles. Within the transmission medium of each paddle is included a transducer array. In a particular embodiment, the array is configured as a monolithic linear assembly. In one embodiment, the array is configured as part of a carriage subsystem for moving the array in an orthogonal direction. Such movement permits the array to scan through the tissue to be examined. In alternative embodiments, the scanning capability may be achieved electronically rather than mechanically. In one embodiment, each array includes 512 transmitting elements and 512 receiving elements. The translation mechanisms of the two acoustic arrays in the upper and lower paddles are configured so that the arrays are positioned substantially opposite each other in the same acoustic plane with the opposing elements facing each other. In a specific embodiment, each of the (say) 1024 transmitting elements sends, in turn, a 2 MHz 3-cycle pulse. The acoustic waveform for each transmitted pulse is received by all (say) 1024 receiving elements. Some of the received signals correspond to transmission through the tissue so that they are received by a receiving element on a different acoustic array than the transmission element. Some of the received signals correspond instead to reflection from the tissue so that they are received by a receiving element on the same acoustic array as the transmission element. The beam pattern produced for each slice is wide in an image plane but narrow in a thickness plane. Once a slice has been taken, both acoustic arrays are moved a distance through the thickness plane, which in one embodiment is approximately 2 mm, and the cycle is repeated. Image reconstruction methods include full aperture tomography, quantitative full aperture tomography, diffraction tomography, full-wave adjoint tomography, full-wave tomography, and a migration algorithm. In one embodiment, the determination of the specific requirements of the migration transformation is made in accordance with the following. The average sound speed in all objects is known to better than about 5% accuracy. The directionality of the sound transmitters and receivers is such that the objects are effectively two-dimensional, and the locus of measurements is a closed curve with the objects located inside the boundary defined by the locus. Thus, in the specific embodiment using the migration transformation, wave traces exist for a fixed set of angle offsets for all positions of the sound transmitter. They form a set of N wave traces for every angle offset along a circular profile. Each of the wave traces permits construction of a position of the reflecting boundary. By increasing the number of wave traces, the accuracy of positioning the reflecting boundary and its diffracting properties are improved.

US 8,532,951 to Roy et al describes an iterative method for calibrating a transducer array characterized by calibration parameters, based on time delay measurements between transducer pairs.

US 8,663,113 to Schmidt et al describes a method for imaging a volume of tissue in a scan region that includes emitting acoustic waveforms toward the volume of tissue, and detecting, with detecting transducers, a plurality of acoustic signals derived from acoustic waveforms interacting with the volume tissue. The method is preferably used to visually characterize tissue of a human breast. The steps of scanning the tissue and detecting acoustic waves are performed during a scan of a patient lying prone on their stomach on a scanner table having a hole opening. The hole in the table allows the breast to protrude through the table and be submerged in an imaging tank filled with water or another suitable fluid as an acoustic coupling medium that propagates acoustic waves. A ring-shaped transducer with transducer elements may be located within the imaging tank and encircle of otherwise surround the breast. Multiple emitters that direct ultrasound pulses toward the tissue, and multiple detectors that receive and record acoustic signals scattering from the tissue, are distributed around the ring transducer. During the scan, the ring transducer passes along the tissue, such as in an anterior-posterior direction between the chest wall and the nipple region of the breast to acquire an acoustic data set, preferably at discrete scanning steps or coronal “slices”. In one preferred embodiment for the ring transducer geometry, an approximation of the Kirchhoff migration algorithm satisfies Dirichlet boundary conditions.

EP 3239706 to Osman et al describes an apparatus with a plurality of transducers which serve as emitter units and receiver units, and allow to emit and receive ultrasound waves respectively. Figure 2 shows an apparatus with a cylindrical measuring area or aperture having a circular base. The object is immersed in a coupling medium and is surrounded by an emitter unit and a receiver unit. Both units are well separated from the measuring object and are supported by a combined emitter/receiver moving unit. The combined emitter/receiver moving unit has the form of a ring and allows two different movements. There is one circular movement of the emitter unit and the receiver unit around the measuring area at one height along the longitudinal axis of the measuring area. The other movement is along the longitudinal axis, i.e. from height to height. At each height along the longitudinal axis various measurements are performed leading to two- dimensional images of the respective slice of the object. The emitter unit is moved by the emitter/receiver moving unit in order to emit the ultrasound waves from various positions. The measurements at the different heights allow to provide a three-dimensional image of the object. Figure 3 shows a measuring area with a cylindrical outer surface with a radius r and a longitudinal axis. The ultrasound transducers serving as receiving units are randomly and/or evenly distributed at the outside of the measuring area. Here, the receiver units are arranged at different heights along the longitudinal axis. The emitter unit is moved by an emitter moving unit along the longitudinal axis to impinge the measuring object within the measuring area with ultrasonic waves from different directions. The receiver units receive ultrasonic waves at various positions after each emission of ultrasonic waves by the emitter unit from different positions relative to the measuring area. Deep knowledge of the ultrasound physics as well as material characteristics (inhomogeneity, anisotropy) is required in order to comprehend what effect will occur within the material. The processed signals are reconstructed to an image of a specified modality. The reflected signals are separately reconstructed using reconstruction methods such as 3-D Synthetic Aperture Focusing Technique and/or iterative reconstruction techniques (Spares Recovery). The transmitted signals are reconstructed in one embodiment via algorithms for through transmission such as Filtered Back Projection, or iterative reconstruction techniques such as Algebraic Iterative Reconstruction Technique.

EP 3906858 to Ruiter et al describes a device for ultrasound-based reflection and transmission tomography. The device comprises a transducer holder configured to hold a plurality of ultrasonic transducers around an imaging volume to be filled with an ultrasonic coupling medium, the transducer holder having an opening for inserting into the imaging volume at least one part of a body to be imaged. The one part of the body to be imaged may be a human breast, particularly a female breast, and images obtained may be subsequently used for diagnosis of breast cancer, for example. The device may provide a “3D ultrasound computer tomography” (3D USCT). 3D may mean an unfocused data-accumulation (e.g. different types of sensors - unfocussed) where the reconstruction can deliver the focusing. The imaging principle of 3D USCT may be based on several thousand ultrasound transducers (e.g. transmitter and receiver) that enclose the breast in an imaging volume filled with an ultrasound coupling medium (e.g. water bath).

Bidikhova and Hall, “Depth Imaging - Above and Beyond Reflections: The Why and How of Your Pre-stack Depth Migration Project,” CSEG Recorder, March 2021, describes the use of pre-stack depth migration (PSDM) in seismic acoustic imaging, and the difficult problems it entails. The article asks:

Why do we even do depth imaging? Would it be practical to ask of someone to spend so much extra time, resources, and money to commission the PSDM work including the associated and complex interval velocity model building? We are told that time migration can be used to solve certain imaging problems, but more complex ones, such as sharp lateral velocity changes, need depth migration. Where do we draw the boundary in the decision between time and depth migration?

■ What are the practical aspects of a PSDM project? Where do we start?

■ How do we gather the data and which data do we have to gather, anyway?

■ Should I involve the interpreters at the onset, or do I need a petrophysicist?

■ What about a geologist?

■ How can I use my well information in PSDM, and do I even need it for all my projects?

■ Can I do without it?

How should I test for migration parameters? ■ How should I set up the migration grid? And what about raytracing?

We will review these questions in the introductory section, and also later on, while discussing the practical aspects of PSDM project execution.

Further questions include, how can I handle anisotropy ambiguity, and what about the near- surface? What are these diffractions, anyways, and why are they so important? And, of course, the never-ending saga of velocity updating: where do I begin, what is my initial velocity model, how do I make it? Which update methods should I choose? Do I need tomography? How can I use my horizons and well tops? What are the rules, anyway, and are there even rules at all? We will try to answer these, and many other questions of practical PSDM project execution. We will also try to touch on selected project management aspects and cover such important topics as quality control and results presentation.

The last and the most important set of questions: Where do we draw the line and pronounce the project and efforts complete and satisfactory? Will it ever be satisfactory? What about interpreters and clients - how can we keep them content and happy? What will happen next, to the data and results, once I pass them on to the interpreter or QI group? How can we feel confident and good about the thorough and intelligent work we have done, while there are so many controversial opinions about the results and the methods themselves?

Additional background art includes US 2006/0287596 to Johnson et al; US 6,005,916 and US 7,684,846 both to Johnson et al; US 2015/0297173 to Klock et al; EP 2148612 to Johnson et al; US 2017/0143304 to Malik et al; US 2016/0210774 to Wiskin et al; US 6,837,854 to Moore et al; US 9,113,835 to Li; US 2013/0041261 to Li et al; US 2013/0204137 to Roy et al; US 2013/0204136 to Duric et al; US 9,763,641 to West et al; US 10,743,837 to Sandhu et al; US 10,368,831 to Duric et al; US 11,071,520 to Sandhu et al; US 11,298,111 to Duric et al; US 2020/0022672 to Duric et al; US 6,786,868 to Stotzka et al; EP 2056124 to Gemmeke et al; EP 2253274 to Gumb et al; EP 2539870 to Stotzka et al; US 8,884,793 to Gemmeke et al; WO 2015/124301 to Gemmeke et al; EP3688486 to Peric et al; EP 4146081 to Zapf; and Vedula et al, “Towards CT-Quality Ultrasound Imaging Using Deep Learning,” arXiv:1710.06304vl [cs.CV] 17Oct2017.

SUMMARY OF THE INVENTION

An aspect of some embodiments of the invention concerns generating a medical ultrasound image of a field of view of a subject’s body, including one or more of 1) using a flexible stretchable array of ultrasound transducers that can wrap around and adapt to bodies of different size and shape, and finding the positions of the transducers before generating the image; 2) using pre-stack depth Kirchhoff migration with transducers situated on multiple sides of the field of view; and 3) taking advantage of the basic similarity of different human bodies to allow automated workflow in taking ultrasound signal data and reconstructing the image, without having to manually intervene.

There is thus provided, in accordance with an exemplary embodiment of the invention, a method of obtaining a medical ultrasound image of a field of view of a subject’s body with a coordinate system comprising grid points, the method comprising: a) placing ultrasound source transducers and ultrasound receiver transducers on the body of the subject, surrounding the field of view; b) recording ultrasound signals each transmitted by one of the source transducers through the body and received by one of the receiver transducer, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) using a velocity model of the field of view to find travel times to each of a plurality of the grid points from each of the source transducers that transmitted signals, and from each of the plurality of grid points to each of the receiver transducers that received signals; d) constructing gathers as a function of grid point and of discrete bins of a parameter related to ray path length, by, for each grid point in the plurality of grid points, and for each of a plurality of bins of the parameter related to ray path length:

1) selecting recorded signals for which an ultrasound ray path going from the source transducer that transmitted the signal to that grid point, and then from that grid point to the receiver transducer that received the signal, has a value of the parameter related to ray path length that is in that bin;

2) evaluating each of the selected recorded signals at its travel time along that ray path, according to the travel times found using the velocity model; and

3) setting the gathers at that grid point and at that bin of the parameter, to a weighted or unweighted sum of the evaluated selected signals; e) correcting the velocity model using information, from dependence of the gathers at a given grid point on value of the parameter, about how ray path length affects travel time; f) correcting the travel times using the corrected velocity model; g) repeating (d), using the corrected travel times; h) iterating (e), (f) and (g) until a convergence criterion is satisfied; i) constructing a pre-stack attribute of the field of view as a function of grid point, using information from dependence of the gathers on grid point and on the parameter; and j) reconstructing an image of the field of view from the pre- stack attribute as a function of grid point.

Optionally, the source transducers and receiver transducers are arranged on the body of the subject such that, for every point in a volume comprising most of the field of view, for every plane that contains the point, there are one or more transducers located on each side of the plane.

Optionally, for most of the source transducers, for every point in a volume comprising most of the field of view, for every plane that contains the point, there are one or more receiver transducers located on each side of the plane, for which an ultrasound signal transmitted by that source transducer and received by that receiver transducer is one of the ultrasound signals recorded.

Optionally, the gathers comprise angle gathers, the parameter related to ray path length is reflection angle at the grid point, and selecting recorded signals comprises finding reflection angle at the grid point for each recorded signal, and selecting signals for which the reflection angle at the grid point is in that bin.

Optionally, finding reflection angle at the grid point for a recorded signal comprises using the velocity model to find a ray path of ultrasound propagating from the source transducer that transmitted the signal to the grid point and from the grid point to the receiver transducer that received the signal.

Optionally, the method comprises using the angle gathers to find pre-stack attributes of one or more of Poisson’s ratio reflectivity, P wave velocity reflectivity, P wave impedance reflectivity, S wave velocity reflectivity, S wave impedance reflectivity, and density reflectivity.

Optionally, constructing the gathers comprises constructing an angle gather for at least one grid point for a reflection angle greater than 60 degrees.

Optionally, correcting the velocity model comprises: a) for the plurality of grid points, finding a value of a trend of the gather with the parameter, for fixed grid point, that would be expected if there were an error in average velocity from an outer surface of the body inward to the grid point; b) using the trend to find a corrected average velocity as a function of grid point; and c) locally inverting the corrected average velocity as a function of grid point to find a corrected interval velocity as a function of grid point.

Optionally, the method also comprises constructing a main reflectivity image by stacking the gathers over values of the parameter.

Optionally, placing the ultrasound source transducers and the ultrasound receiver transducers on the body of the subject comprises wrapping around a part of the body a flexible stretchable frame, adapted to conform to different body shapes, with the transducers mounted on the frame and held in acoustic contact with the body by the frame.

Optionally, the method also comprises obtaining information, about relative positions, at times when the ultrasound signals are transmitted or received, of the source transducers and the receiver transducers.

Optionally, the coordinate system comprises an (r, θ, z) coordinate grid, and constructing gathers comprises constructing the gathers with a discrete set of grid coordinates (θ, z) defining the locations of gathers extending into the field of view away from an outer surface of the subject’s body, and wherein there is a range of values of θ and a range of values of z that include at least ten pairs (θ, z) that belong to the discrete set, and a range of values of r such that all of the gathers defined by said pairs (θ, z) extend over said range of r, with r increasing along each of said gathers in a direction going into the field of view further away from a closest outer surface of the subject’s body, and wherein, over said range of r, the average distance between said gathers at a given r is smaller for larger values of r.

Optionally, the field of view has a tube-shaped outer surface, corresponding to an outer surface of the subject’s body, and wherein the coordinate system comprises an (r, θ, z) coordinate grid, with θ a periodic azimuthal coordinate around the tube, z an axial coordinate along the tube, and r a coordinate going into the field of view from the outer surface, and wherein constructing the gathers comprises constructing the gathers with a discrete set of grid coordinates (θ, z) defining the locations of the gathers extending into the field of view away from the outer surface, and wherein there is a range of values of z that include at least ten pairs (θ, z) that belong to the discrete set, and a range of values of r such that all of the gathers defined by said pairs (θ, z) extend over said range of r, with r increasing along each of said gathers in a direction going into the field of view away from the outer surface, and wherein, for said range of z and said range of r, the field of view extends over a full range of values of 0.

Optionally, correcting the velocity model using the gathers comprises: a) finding a trend in a dependence of the gather on the parameter for fixed r, θ and z, that would be expected if there were an error in the velocity model, and using the trend to find an average velocity as a function of r for each pair (θ, z) that corrects the error; and b) finding a corrected interval velocity as a function of r, θ, and z, at least over said ranges of r, θ, and z, by locally inverting the average velocity as a function of r for each pair (θ, z), and smoothing the interval velocity, as a function of r and z, at least over 0. Optionally, the parameter comprises offset between source transducer and receiver transducer.

Optionally, the method also comprises using the pre-stack attribute to find a corresponding post-stack attribute of at least a portion of the field of view as a function of grid point.

Optionally, the pre-stack attribute comprises a reflectivity of a quantity, the post-stack attribute corresponding to that pre- stack attribute is the quantity, and finding the corresponding post-stack attribute comprises performing a local inversion procedure using an optimization algorithm.

Optionally, using the pre-stack attribute to find a corresponding post-stack attribute comprises finding one or more of Poisson’s ratio, P wave velocity, P wave impedance, S wave velocity, S wave impedance, and density.

Optionally, the method comprises finding one or more post-stack attributes by performing a mathematical manipulation of a pre-stack attribute or of a post-stack attribute.

Optionally, performing the mathematical manipulation comprises taking a complex trace.

There is further provided, according to an exemplary embodiment of the invention, a method of obtaining at least a first medical ultrasound image of a field of view of a first body part of a first subject, the method comprising: a) wrapping a flexible stretchable frame, with ultrasound source transducers and ultrasound receiver transducers mounted on it, around the first body part, with the ultrasound transducers held by the frame in good acoustic contact with an outer surface of the first body part; b) recording a first set of ultrasound signals, each signal transmitted from one of the source transducers through the first body part and received by one of the receiver transducers, while the transducers are held in good acoustic contact with the outer surface of the first body part, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) determining, from the recorded first set of ultrasound signals, information about relative spatial positions of the source transducers that transmitted one of the first set of signals and the receiver transducers that received one of the first set of signals, at times when the first set of ultrasound signals were transmitted or received; and d) reconstructing the first ultrasound image from the recorded first set of ultrasound signals and the information about the relative spatial positions of the transducers at times when the signals were transmitted or received, after determining the information about the relative spatial positions of the transducers. Optionally, reconstructing the first ultrasound image from the recorded first set of ultrasound signals uses a coordinate system comprising a plurality of grid points, and comprises: a) using a velocity model of the field of view to find travel times to each grid point from each source transducer that transmitted a signal and from each grid point to each receiver transducer that received a signal; and b) constructing gathers of the recorded signals, using the travel times, as a function at least of grid point; and wherein constructing the gathers is done after determining the information about the relative spatial positions of the transducers.

Optionally, the method obtains also a second medical ultrasound image of a field of view of a second body part, of the first subject or a different subject, the second body part having a shape, size, or both that differs from the first body part, the method also comprising: a) after recording the first set of ultrasound signals, wrapping the frame around the second body part, with the ultrasound transducers held by the frame in good acoustic contact with an outer surface of the second body part, adjusting the frame if necessary to produce the good acoustic contact; b) recording a second set of ultrasound signals each signal transmitted from one of the source transducers through the second body part and received by one of the receiver transducers, while the transducers are held in good acoustic contact with the outer surface of the second body part, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) determining, from the recorded second set of ultrasound signals, information about relative spatial positions of the source transducers that transmitted one of the second set of signals and the receiver transducers that received one of the second set of signals, at times when the second set of ultrasound signals were transmitted or received; and d) calculating the second ultrasound image from the recorded second set of ultrasound signals and the information about the relative spatial positions of the transducers at times when the signals were transmitted or received, after determining the information about the relative spatial positions of the transducers.

Optionally, the first body part is a trunk of the first subject, and the transducers are mounted on the frame, and the frame is wrapped around the trunk, such that, along at least one path along the frame encircling the trunk, there is no gap between transducers that is wider than 20% of the length of the path. Optionally, the field of view comprises at least 50% of at least one axial cross-section of the trunk.

Optionally, the field of view is three dimensional, and includes at least 50% of the volume of one or more of: 1) the lungs and the heart; 2) the stomach and the liver; 3) the small intestine and the colon.

Optionally, the method obtains also a second medical ultrasound image of a field of view of a trunk of a second subject, differing in size, shape, or both from the trunk of the first subject, the method comprising: a) after recording the first set of ultrasound signals, wrapping the frame around the trunk of the second subject, with the ultrasound transducers held by the frame in good acoustic contact with an outer surface of the trunk of the second subject, adjusting the frame if necessary to produce the good acoustic contact; b) recording a second set of ultrasound signals, each signal transmitted from one of the source transducers through the trunk of the second subject and received by one of the receiver transducers, while the transducers are held in good acoustic contact with the outer surface of the trunk of the second subject, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) determining, from the recorded second set of ultrasound signals, information about relative spatial positions of the source transducers that transmitted one of the second set of signals and the plurality of receiver transducers that received one of the second set of transducers, at times when the second set of ultrasound signals were transmitted or received; and d) calculating the second ultrasound image from the recorded second set of ultrasound signals and the information about the relative spatial positions of the transducers at times when the signals were transmitted or received, after determining the information about relative spatial positions of the transducers.

Optionally, determining information about relative spatial positions of the transducers comprises determining spatial positions of the transducers with a standard deviation error of less than 0.5 mm.

Optionally, the method also comprises providing initial estimates for the relative spatial positions of the source transducers that transmitted one of the signals in the first set of signals and the receiver transducers that received one of the signals in the first set of signals, at the times when the first set of ultrasound signals were transmitted or received, wherein determining information about the relative spatial positions of the transducers comprises using the initial estimates. Optionally, providing initial estimates for the relative spatial positions of the transducers comprises: a) providing information on locations where the transducers are mounted on the frame; and b) using the information on the locations where the transducers are mounted on the frame, and a model of the first subject’s body, to estimate relative spatial positions of the transducers when the frame is wrapped around the first body part.

Optionally, the frame comprises one or more sensors that generate data of a distance between the transducers in one or more pairs of transducers mounted adjacent to each other on the frame, one or more sensors that generate data of a curvature or angle between the transducers in one or more sets of three transducers mounted adjacent to each other on the frame, or both.

Optionally, providing the initial estimates for the relative spatial positions of the transducers comprises using the data generated by said sensors.

Optionally, providing the initial estimates for the relative spatial positions of the transducers also comprises: a) providing information on locations where the transducers are mounted on the frame; and b) using the information on the locations where the transducers are mounted on the frame, and a model of the first subject’s body, to estimate relative spatial positions of the transducers when the frame is wrapped around the first body part.

Optionally, determining information about the relative spatial positions of the transducers from the recorded first set of ultrasound signals comprises determining information about the relative spatial positions of the transducers from only a sub-set of the recorded first set of ultrasound signals where the source transducer and the receiver transducer are located relatively close to each other, compared to distances between the source transducer and the receiver transducer for ultrasound signals in the recorded first set that are not in the sub-set.

Optionally, determining information about relative spatial positions of the transducers comprises using a computer running a program configured to improve its performance by machine learning.

Optionally, determining information about relative spatial positions of the transducers comprises using a computer running a neural network configured to determine the positions from the recorded first set of ultrasound signals.

Optionally, the neural network comprises a convolutional neural network (CNN).

Optionally, the CNN has a U-shaped architecture. Optionally, running the CNN comprises: a) receiving a primary input data volume comprising the recorded signals from the first set of source transducer-receiver transducer pairs; b) passing the primary input data volume through the U-shaped architecture to generate a first intermediate data volume; c) receiving a secondary input data volume comprising initial estimates of the positions of the source transducers and the receiver transducers in the first set of pairs; d) performing a first combining of data volumes using the first intermediate data volume and the secondary input data volume, by one of:

1) element-by-element addition of the two data volumes, or element-by-element addition of the two data volumes after first passing one or both of them at least through a dense layer;

2) element-by-element multiplication of the two data volumes, or element=by- element multiplication of the two data volumes after first passing one or both of them at least through a dense layer; and

3) concatenation of the two data volumes; to generate a second intermediate data volume; e) using the second intermediate data volume to generate an output data volume comprising positions of the source transducers and the receiver transducers in the first set of pairs.

Optionally, using the second intermediate data volume to generate the output data volume comprises performing a second combining of data volumes using the second intermediate data volume and one of the first intermediate data volume and the secondary input data volume, by one of:

3) concatenation of the two data volumes; to generate a third intermediate data volume, and using the third intermediate data volume to generate the output data volume.

Optionally, using the third intermediate data volume to generate the output data volume comprises performing a third combining of data volumes using the third intermediate data volume and a different one of the first intermediate data volume and the secondary input data volumes than used in the second combining of data volumes, by one of:

3) concatenation of the two data volumes; to generate a fourth intermediate data volume, and using the fourth intermediate data volume to generate the output data volume.

Optionally, the first combining of data volumes, the second combining of data volumes, and the third combining of data volumes each use a different one of element-by-element addition, element-by-element multiplication, and concatenation.

Optionally, each of the source transducers transmits signals to the receiver transducers at a different time, and the output data volume comprises a position of each of the receiver transducers at each of the different times when one of the source transducers transmits signals to it.

There is further provided, according to an exemplary embodiment of the invention, a method of training a neural network for use in an exemplary embodiment of the invention, or for training a program for use in an exemplary embodiment of the invention, to determine the relative spatial positions of the transducers from the recorded first set of ultrasound signals, the method comprising using numerical solutions of a 3-D wave equation for ultrasound waves transmitted through a body part from a source transducer to a receiver ultrasound transducer to simulate the recorded first set of ultrasound signals, using at least one model of ultrasound propagation properties as a function of position in the body part, for a plurality of different combinations of positions of source transducers and positions of receiver transducers on the outer surface of the body part.

Optionally, using at least one model of ultrasound propagation properties as a function of position in the body part comprises using a plurality of models of ultrasound propagation properties as a function of position in the body part and shape and size of the body part, corresponding to a range of subjects in a population of subjects.

There is further provided, in accordance with an exemplary embodiment of the invention, a method of obtaining one or more types of medical ultrasound images of bodies of subjects belonging to a population, the method comprising: a) providing a subject belonging to the population, general information about the subject, and a type of ultrasound image to be obtained of a field of view of the subject’s body; b) placing on the subject’s body an array comprising ultrasound source transducers and ultrasound receiver transducers, the transducers in good acoustic contact with an outer surface of the subject’s body, and positioned so that the receiver transducers can receive ultrasound signals transmitted into the field of view by the source transducers; c) automatically providing an initial velocity model for sound speed as a function of position in at least a part of the subject’s body including the field of view, one or more parameters for transmitting a set of ultrasound signals through the subject’s body from a plurality of the source transducers to a plurality of the receiver transducers, one or more parameters for finding travel times using a velocity model, one or more parameters for constructing gathers from recorded received ultrasound signals, one or more parameters for correcting the velocity model, one or more parameters for iterating the correction of the velocity model, one or more parameters for constructing a pre-stack attribute from the gathers, and one or more parameters for reconstructing an ultrasound image from the pre-stack attribute; d) recording a set of ultrasound signals, each signal transmitted by one of the source transducers through the subject’s body and received by one of the receiver transducers, according to the parameters for transmitting the set of ultrasound signals through the subject’s body, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; e) using a velocity model of the field of view to find travel times to each of a plurality of grid points from each source transducer that transmitted one of the signals and from each of the plurality of grid points to each receiver transducer that received one of the signals, according to the parameters for using a velocity model to find travel times; f) constructing gathers as a function of grid point and of discrete bins of a parameter related to ray path length, by, for each grid point in the plurality of grid points, and for each of a plurality of bins of the parameter related to ray path length:

2) evaluating each of the selected recorded signals at its travel time along that ray path, according to the travel times found using the velocity model; and 3) setting the gathers at that grid point and at that bin of the parameter, to a weighted or unweighted sum of the evaluated selected signals; according to the parameters for constructing gathers; g) correcting the velocity model using information, from dependence of the gathers at a given grid point on value of the parameter, about how ray path length affects travel time, according to the parameters for correcting the velocity model; h) correcting the travel times using the corrected velocity model, according to the parameters for using a velocity model to find travel times; i) repeating (f), using the corrected travel times; j) iterating (g), (h) and (i) until a convergence criterion is satisfied, according to the parameters for iterating the correction of the velocity model and the gathers; k) constructing a pre-stack attribute of the field of view as a function of grid point, using information from dependence of the gathers on grid point and on the parameter, according to the parameters for constructing a pre-stack attribute; and l) reconstructing an image of the field of view from the pre- stack attribute as a function of grid point, according to the parameters for reconstructing an ultrasound image from the pre- stack attribute.

Optionally, the array comprises a flexible stretchable frame on which the transducers are mounted, and placing the array on the body comprises placing the flexible stretchable frame around the body.

Optionally, one or more of automatically providing the initial velocity model, automatically providing one or more parameters for transmitting a set of ultrasound signals through the subject’s body from a plurality of the source transducers to a plurality of the receiver transducers, one or more parameters for finding travel times using a velocity model, one or more parameters for constructing gathers from recorded received ultrasound signals, one or more parameters for correcting the velocity model, one or more paramaters for iterating the correction of the velocity model, one or more parameters for constructing a pre-stack attribute from the gathers, and one or more parameters for reconstructing an ultrasound image from the pre- stack attribute, comprises automatically providing based on one or more of general information about the subject, the type of ultrasound image, and data about spacing of the transducers obtained automatically after placing the frame around the subject’s body.

Optionally, the method also comprises automatically providing one or more parameters for finding reflection angle, wherein the parameter related to ray path length comprises reflection angle, and selecting recorded signals comprises finding reflection angle at the grid point for each recorded signal, and selecting signals for which the reflection angle at the grid point is in that bin.

Optionally, the method also comprises automatically providing one or more parameters for using a velocity model to find ray paths of propagating ultrasound, where finding reflection angle, at the grid point, of an ultrasound signal travelling from a source transducer to the grid point to a receiver transducer comprises using the velocity model to find a ray path of ultrasound propagating from the source transducer to the grid point and from the grid point to the receiver transducer, according to the parameters for using the velocity model to find ray paths.

Optionally, constructing gathers comprises constructing gathers on an (r, θ, z) coordinate grid, with a discrete set of grid coordinates (θ, z) defining the locations of gathers extending into the field of view away from an outer surface of the subject’s body, and wherein there is a range of values of θ and a range of values of z that include at least ten pairs (θ, z) that belong to the discrete set, and a range of values of r such that all of the angle gathers defined by said pairs (θ, z) extend over said range of r, with r increasing along each of said angle gathers in a direction going into the field of view away from a closest outer surface of the subject’s body, and wherein, over said range of r, the average distance between said angle gathers at a given r is smaller for larger values of r.

Optionally, correcting the velocity model using the gathers comprises: a) finding a trend in a dependence of the gather on the parameter related to path length, for fixed r, θ and z, that would be expected if there were an error in the velocity model, and using the trend to find an average velocity as a function of r for each pair (θ, z) that corrects the error; and b) finding a corrected interval velocity as a function of r, θ, and z, at least over said ranges of r, θ, and z, by locally inverting the average velocity as a function of r for each pair (θ, z), and smoothing the interval velocity, as a function of r and z, at least over 0.

There is further provided, according to an exemplary embodiment of the invention, a system for obtaining a medical ultrasound image of a field of view of a body part of a subject with a coordinate grid comprising grid points, the system comprising: a) an ultrasound array comprising:

1) a flexible stretchable frame of a two-dimensional structure;

2) a plurality of ultrasound source transducers mounted on the frame in a distributed arrangement over an area of the two-dimensional structure;

3) a plurality of ultrasound receiver transducers, one or more of which are also source ultrasound transducers or none of which are also source ultrasound transducers, mounted on the frame in a distributed arrangement over an area of the two-dimensional structure; and 4) a fastener to hold different ends of the two-dimensional structure of the frame together, holding the frame in place when the frame is wrapped around the body part; wherein the transducers are mounted on the frame with faces facing a same direction normal to the area of the two-dimensional structure when the frame is unwrapped and flat, such that the transducers are held pressed in good acoustic contact with an outer surface of the body part when the frame is wrapped and fastened around the body part in at least one configuration, the array adapted for wrapping and fastening around the body part of the subject with the transducers held pressed in good acoustic contact with an outer surface of the body part; b) a controller, comprising an output interface, configured to:

1) cause each of a plurality of the source transducers to transmit an ultrasound signal into the body part at different times;

2) cause one or more of the receiver transducers to record the ultrasound signals received from each of the plurality of source receivers, and to store data of each of the recorded signals in association with an identifier of its source transducer and an identifier of its receiver transducer;

3) use a velocity model of the field of view to find travel times to each of a plurality of the grid points from each of the source transducers that transmits signals and from each of the plurality of grid point to each receiver transducer that records signals;

4) construct gathers as a function of grid point and of discrete bins of a parameter related to ray path length, by, for each grid point in the plurality of grid points, and for each of a plurality of bins of the parameter related to ray path length: i) selecting recorded signals for which an ultrasound ray path going from the source transducer that transmitted the signal to that grid point, and then from that grid point to the receiver transducer that received the signal, has a value of the parameter related to ray path length that is in that bin; ii) evaluating each of the selected recorded signals at its travel time along that ray path, according to the travel times found using the velocity model; and iii) setting the gathers at that grid point and at that bin of the parameter, to a weighted or unweighted sum of the evaluated selected signals;

5) correct the velocity model using information, from dependence of the gathers at a given grid point on value of the parameter, about how ray path length affects travel time;

6) correct the travel times using the corrected velocity model;

7) repeat (4), using the corrected travel times;

8) iterate (5), (6) and (7) until a convergence criterion is satisfied; 9) construct a pre-stack attribute of the field of view as a function of grid point, using information from dependence of the gathers on grid point and on the parameter;

10) reconstruct an image file of the field of view from the pre- stack attribute as a function of grid point; and

11) output the ultrasound image file through the output interface; and c) one or more power supplies that power the source transducers, the receiver transducers, and the controller.

Optionally, the ultrasound array also comprises sensors incorporated in the frame that measure distance along the frame between one or more pairs of the transducers, angles in the two- dimensional structure of the frame between one or more sets of three of the transducers, or both, when the frame is wrapped around and the body part and fastened with the transducers held pressed against the outer surface of the body part.

There is further provided, in accordance with an exemplary embodiment of the invention, a method of obtaining a medical ultrasound image of a field of view of a subject’s body with grid points, comprising: a) placing ultrasound source and receiver transducers on the body, surrounding the field of view; b) recording ultrasound signals transmitted by a source transducer through the body and received by a receiver transducer; c) using a velocity model to find travel times between grid points and source and receiver transducers; d) constructing gathers as a function of grid point and a parameter related to ray path length, using pre-stack depth Kirchhoff migration; e) correcting the velocity model using information from the gathers about how ray path length affects travel time; f) constructing a pre- stack attribute of the field of view as a function of grid point, from the gathers; and g) reconstructing an image of the field of view from the pre- stack attribute as a function of grid point.

Unless otherwise defined, all technical and/or scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the invention pertains. Although methods and materials similar or equivalent to those described herein can be used in the practice or testing of embodiments of the invention, exemplary methods and/or materials are described below. In case of conflict, the patent specification, including definitions, will control. In addition, the materials, methods, and examples are illustrative only and are not intended to be necessarily limiting.

Implementation of the method and/or system of embodiments of the invention can involve performing or completing selected tasks manually, automatically, or a combination thereof. Moreover, according to actual instrumentation and equipment of embodiments of the method and/or system of the invention, several selected tasks could be implemented by hardware, by software or by firmware or by a combination thereof using an operating system.

For example, hardware for performing selected tasks according to embodiments of the invention could be implemented as a chip or a circuit. As software, selected tasks according to embodiments of the invention could be implemented as a plurality of software instructions being executed by a computer using any suitable operating system. In an exemplary embodiment of the invention, one or more tasks according to exemplary embodiments of method and/or system as described herein are performed by a data processor, such as a computing platform for executing a plurality of instructions. Optionally, the data processor includes a volatile memory for storing instructions and/or data and/or a non-volatile storage, for example, a magnetic hard-disk and/or removable media, for storing instructions and/or data. Optionally, a network connection is provided as well. A display and/or a user input device such as a keyboard or mouse are optionally provided as well.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING(S)

Some embodiments of the invention are herein described, by way of example only, with reference to the accompanying drawings. With specific reference now to the drawings in detail, it is stressed that the particulars shown are by way of example and for purposes of illustrative discussion of embodiments of the invention. In this regard, the description taken with the drawings makes apparent to those skilled in the art how embodiments of the invention may be practiced.

In the drawings:

FIG. 1 schematically shows a perspective view of an array of acoustic sources and geophones on the surface of the Earth, used for seismic imaging of a field of view beneath the surface, according to the prior art;

FIG. 2 schematically shows a system, including an ultrasound transducer array, for medical ultrasound imaging, according to an exemplary embodiment of the invention; FIG. 3 is a flowchart for a method of generating a medical ultrasound image using the system of FIG. 2, for example, according to an exemplary embodiment of the invention;

FIG. 4 is a flowchart showing further details of the method of FIG. 3, according to an exemplary embodiment of the invention;

FIGS. 5 A and 5B are schematic cut-away perspective views showing a cylindrical body part with source and receiver transducers mounted at different positions on the body surface, showing sources that transmit ultrasound to an imaging location inside the body and receivers that detect the ultrasound that reflects from the imaging location, with different reflection angles and different imaging directions;

FIG. 6 is a schematic perspective view of an elliptical cylindrical coordinate system (r, θ, z), for a field of view of a body part with a surface that is roughly an elliptical cylinder, used to generate an ultrasound image of the field of view, for example with the method of FIGS. 3 and 4, according to an exemplary embodiment of the invention;

FIG. 7A is a flowchart showing further details of the method of FIGS. 3 and 4, according to an exemplary embodiment of the invention;

FIG. 7B schematically shows an image point in a field of view used for seismic imaging, with an acoustic source and receiver at positions for which the reflection angle would be 80°, and an image point in a field of view used for ultrasound medical imaging, with a source and receiver at positions for which the reflection angle would be 80°;

FIGS. 8 A and 8B schematically show perspective views of an ultrasound array similar to that shown in FIG. 2, showing two different local sets of nearby source and receiver transducers used for finding precise positions of the transducers for example according to the method of FIG. 7A, using a program configured for being trained by machine learning,;

FIG. 9A is a block diagram showing the architecture of a convolutional neural network that can be trained by machine learning for finding precise positions of the transducers, for example for use in the method of FIG. 7A; and

FIG. 9B is a block diagram showing further details of one module in FIG. 9A;

FIG. 10A is a block diagram showing how training examples are generated for training a convolutional neural network, such as the one diagrammed in FIG. 9A, used for finding precise positions of transducers, for example according to the method of FIG. 7A, according to an exemplary embodiment of the invention; and

FIG. 10B is a flowchart for a method of training such a convolutional neural network, using training examples such as those described in FIG. 10A; FIG. 11A is a schematic cross-sectional side view showing source and receiver ultrasound transducers mounted on a surface of a body part being imaged, for example using the system of FIG. 2 and the method of FIGS. 3, 4, and 7A, showing ray paths followed by ultrasound signals reflecting with different reflection angles that contribute to the migrated reflectivity from a same location in an angle gather in an idealized field of view with a flat surface and a uniform sound velocity;

FIG. 1 IB schematically shows the data in an angle gather as a function of travel time t and reflection angle a, for a field of view with only a single point feature that has high reflectivity at a specific depth, for example obtained using the system of FIG. 2 and the method of FIGS. 3, 4 and 7A, using a correct velocity model and using an incorrect velocity model, according to an exemplary embodiment of the invention;

FIG. 12A schematically shows data in an angle gather, more realistic than the data shown in FIG. 1 IB, showing reflectivity amplitude M as a function of vertical travel time / and reflection angle a, for example obtained using the system of FIG. 2 and the method of FIGS. 3, 4 and 7A, according to an exemplary embodiment of the invention;

FIG. 12B schematically shows the reflectivity amplitude M as a function of reflection angle a, for vertical travel time t = 1110 ms, according to the data shown in FIG. 12A, showing how M(t, a) is fitted to a form I + Gsin²a for each value of / in the angle gather, to find a correction to the velocity model, for example according to the method of FIG. 7A, according to an exemplary embodiment of the invention; and

FIG. 13 is a flowchart showing further details of calculating pre-stack attributes and poststack attributes for a medical ultrasound image, according to the method of FIG. 4.

DESCRIPTION OF SPECIFIC EMBODIMENTS OF THE INVENTION

The present invention, in some embodiments thereof, relates to an ultrasound imaging system and method, and, more particularly, but not exclusively, to a medical ultrasound system and method for generating an ultrasound image of a field of view of the interior of a human body.

The methods of ultrasound medical imaging described herein grew out of methods used for geological acoustic imaging in the oil, gas, and mineral industries. As used herein, “geological acoustic imaging,” “geological imaging,” “seismic acoustic imaging,” and “seismic imaging” have the same meaning. People have invented several ways to reconstruct acoustic images of geological formations, including post- stack migration, pre- stack migration, including Kirchhoff migration and reverse time migration, and FWI (Full Wave Inversion). Migration is a process by which reflecting events are re-located in either space or time to the location where the reflecting event had occurred rather than the location where it was recorded at the surface, thereby creating a more accurate and focused image, Focusing is achieved because all scattered energy returns to the point it was scattered from. In reverse time migration, waves are followed backward in time to the point they scattered from by solving the wave equation. In Kirchhoff migration, recorded acoustic signals, from different combinations of acoustic sources and acoustic receivers, are added up, each evaluated at a time that corresponds to the travel time of a sound wave travelling from that source to that receiver, and reflecting from a point (x, y, z) in the field of view, and optionally with each signal weighted by a scale factor. Although each signal also includes contributions due to reflections from other points in the field of view, those contributions tend to cancel out due to destructive interference, while the contributions due to reflection from that point (x, y, z) tend to reinforce by constructive interference, so the sum of the time-delayed signals provides information about acoustic reflectivity at that point (x. y, z). Kirchhoff migration can be post-stack Kirchhoff migration, which may provide information only about reflectivity, or pre-stack Kirchhoff migration, in which the summed signals are sorted by some parameter, such as reflection angle at the point (x, y, z), which can be used to provide information also about other acoustic properties at (x, y, z), such as sound speed, shear wave speed, and Poisson’s ratio. FWI is a process that aims to find a high resolution velocity model from the recorded data. Unlike migration, FWI does not attempt to provide a reflectivity image, but rather a velocity image. FWI solves a large inversion matrix that contains the differences between the recorded data and simulated data (using direct solution of the wave equation as a function of the velocity model), and finds the best velocity model that minimizes this difference. Further details about migration and FWI are provided in some of the references cited above in the “Field and Background” section, and in the rest of the “Description of Specific Embodiments” below.

As computers improved, more sophisticated methods became possible. But given the value of large deposits of oil and other minerals, it was possible to use methods that were very labor intensive and required a lot of manual decision-making each time a new geological formation was imaged. As described in the paper by Bidikhova and Hall cited above, it is well known among experts in geological acoustic imaging that pre-stack depth Kirchhoff migration with velocity model building is very labor intensive. In particular, in this method recorded acoustic signals can be added together into “angle gathers” M(x, y, z, a), where (x, y, z) represents the spatial coordinates of a point in the underground field of view, and a represents the angle of reflection of a sound wave (in the ray optics approximation) reflecting from that point. If an accurate velocity (sound speed) model is being used to calculate the ray paths of the sound waves in the field of view, then we do not expect to see a systematic trend of M(x, y, z, a) with a, for fixed x, y and z. If a trend is observed, for some values of (x, y, z), then it is assumed to be due to an error in the velocity model, and it can be used to correct the velocity model over the field of view, which in turn can be used to find corrected angle gathers M(r, θ, z, a). This process is then iterated until it converges on a correct velocity model and reflectivity image of the field of view. It is this process that is very labor intensive in geological imaging, requiring extensive manual decision making that depends on the details of the geological formation and on the locations and characteristics of the acoustic sources and sensors, and on the protocol for acquiring data. The cost of imaging one geological formation could be as high as a million dollars, but it is still worthwhile in the oil, gas and mining industry. Even if the manual decisions that have to be made are not very timeconsuming for an expert to make, just the fact that such manual decisions have to be made at all means that the workflow cannot be automated, and the method may not be practical for medical imaging.

When the idea developed to adapt geological acoustic imaging methods to ultrasound medical imaging, the people who developed the methods were necessarily experts in geological acoustic imaging. Because those experts knew that pre-stack depth Kirchhoff migration, including building the velocity model, was very computationally expensive as well as labor intensive, and that medical ultrasound images should each cost no more than hundreds of dollars to produce including labor costs, they did not seriously consider using pre-stack depth Kirchhoff migration, with velocity model correction, for ultrasound medical images as a production tool. The most commonly used geology-inspired method for medical ultrasound imaging has been FWI, and poststack migration, as well as pre-stack time migration methods, are also sometimes used. FWI is a computationally expensive method, but the velocity analysis is implicit in it, rather than being a separate step, so it does not require repeated manual decision making.

It was an inventive insight of the present inventor that in medical ultrasound imaging, prestack depth Kirchhoff migration workflow, including correcting the velocity model, can be fully automated, greatly bringing down the cost and making it economically feasible. This is possible because people’s bodies are much more similar to each other than geological formations, so it is possible to design a single acquisition device and a single protocol that can be used to image the same body part of many different patients.

Once this is understood, it becomes possible to take advantage of the many desirable features that pre-stack depth Kirchhoff migration has for ultrasound medical imaging, both in comparison with pre-stack depth Kirchhoff migration for geological acoustic imaging, and in comparison with ultrasound medical imaging using FWI and post-stock migration methods. These advantages may be particularly important for imaging a field of view with large volume and high resolution, and using a large array of closely spaced transducers. For such ultrasound medical images, pre-stack depth Kirchhoff migration may not only be practical, but may be by far the best method to use.

In geological acoustic imaging, the field of view is a large volume located underground, down to a large depth, but the transducers for emitting and receiving sound waves are all on or very close to the surface of the earth. So the transducers are located on only one side of the field of view, on its top surface. In medical ultrasound imaging, transducers can be located on multiple sides of the field of view, optionally completely surrounding it, or surrounding it in all directions azimuthally, for example surrounding the trunk of the body if the field of view is all or part of the trunk of the body. For pre-stack depth Kirchhoff migration, this has several advantages. In geological imaging, where the area on the surface of the earth over which transducers are distributed is typically not much wider than the depth of the field of view, reflection angles are limited, typically to no more than 40 degrees. In ultrasound medical imaging, where the body can be completely surrounded by transducers, reflection angles can extend up to 90 degrees, or close to 90 degrees. This is important with pre- stack depth Kirchhoff migration, where sorting data by reflection angle is used for determining the velocity model (sound speed as a function of position in the field of view), as well as for measuring other attributes as a function of position in the field of view, such as Poisson ratio, elasticity, and density, all of which can be important for distinguishing different types of tissue, including healthy and diseased tissue, in medical imaging.

In geological acoustic imaging, including with pre-stack depth Kirchhoff migration, there are edge effects in the part of the field of view near its lateral edges, which can degrade the image over a large part of the field of view. In medical ultrasound imaging, if transducers are located 360 degrees around the trunk, or around a limb, azimuthally, then there will be no edge effects, at least azimuthally, and more of the field of view will have a good quality image, though there might still be edge effects in an axial direction.

In geological acoustic imaging, using pre-stack depth Kirchhoff migration, angle gathers generally extend vertically, in the z direction with Cartesian coordinates, perpendicular to the surface of the earth, down into the field of view. The angle gathers are all parallel to each other, and signal data from deeper down in the field of view is weaker and noisier, and provides less accurate information about the sound velocity at those locations. This is also true because at greater depths beneath the surface, the usable range of reflection angles is smaller, and determining sound velocity depends on measuring how signals vary with reflection angle. As a result, the image at greater depths tends to be less accurate, and especially the sound velocity tends to be less accurate. In medical ultrasound imaging, if there are transducers located over 360 degrees, or over any wide range of angles azimuthally, for example around the trunk, angle gathers will still extend into the field of view perpendicular to the surface, in this case in the r direction with cylindrical coordinates. Although the signal data for any given angle gather will be weaker and noisier deeper into the body, the angle gathers will get closer together, the deeper they go into the body. As a result, the image attributes in neighboring angle gathers will resemble each other more, the deeper they go into the body, and neighboring angle gathers can be used to correct each other. This can compensate for the weaker signal and greater noise within each angle gather, at greater distances into the body, and the image quality further into the body need not be worse than near the surface. And the usable range of reflection angles will be no less deep inside the body than near the surface, so determination of the sound velocity will not tend to be worse deep inside the body.

Another potential advantage of having transducers surrounding the field of view in a cylindrical body, especially the trunk, is that, for any part of the surface of any organ, there will be ultrasound waves propagating approximately normal to it, from nearby source transducers, and reflecting from the surface of the organ toward nearby receiver transducers, producing strong signals that can be used to reconstruct good images of that part of the surface of the organ. In seismic imaging, on the other hand, well below the ground sound waves will be propagating in only a limited range of directions around the vertical. For geological structures that have approximately vertical surfaces, those sound waves will reflect from those surfaces at a glancing angle, for which the reflectivity may be low, so the signals may be weak and the resulting images may be poor quality.

Yet another potential advantage of having transducers surrounding the field of view in a cylindrical body part is that the acoustic power per area illuminating a given point in the field of view, coming from transducers in multiple directions, will be greater, relative to the power per area at the transducer, than if transducers are located on only one side of the field of view, as in geological acoustic imaging. This can result in a better quality image, or a faster image acquisition time, for a given power per area of the transducers.

Using pre-stack depth Kirchhoff migration for medical ultrasound imaging also has important advantages over using FWI for medical ultrasound imaging, for a given field of view and resolution. With pre-stack depth Kirchhoff migration, the calculation of the image is broken into a set of independent steps, one following another, and each calculation may be simpler than with FWI. First, the receiver data is combined into gathers, such as angle gathers. Then the angle gathers are used to correct the velocity model. Then the corrected gathers are stacked to form a high resolution reflectivity image. The gathers can then be used to form high resolution images of other attributes, such as sound speed, at much higher resolution than the velocity model, as well as Poisson’s ratio, density, and elastic modulus. When FWI is used, a single large inversion calculation is made, for example inverting a large matrix, to find everything in the image, in particular to find a high resolution map of the sound velocity. Especially for a large field of view at high resolution, where there are a very large number of grid points at which the image is being found, and a very large number of transducers, the calculation of the sound velocity at each grid point may tend to be noisy and unstable, finding local minima in the quantity that is being optimized, yielding an inaccurate image, rather than finding the global minimum that is being sought. It should be noted that the velocity model calculated in pre- stack depth Kirchhoff migration, which is needed for finding the ray optics paths and travel times of the sound waves, can be very low resolution, even for a high resolution reflectivity image. With pre-stack depth Kirchhoff migration, a high resolution image of the sound speed may be calculated later in a separate step, after the high resolution reflectivity image has been found. That calculation will be a local inversion calculation for each grid point, rather than a global inversion calculation as in FWI, so will be orders of magnitude cheaper computationally, and will not suffer from instability.

In geological acoustic imaging, acoustic sources and receivers are located in known fixed positions on or near the surface of the earth. Similarly, in many applications of geological acoustic methods to medical ultrasound imaging, the ultrasound transducers are located in a rigid frame at fixed positions, for example mounted on the head for imaging the brain, or mounted on the outside of a rigid fluid-filled container for imaging a breast, or other soft body parts. Having the transducers at known fixed positions, relative to the field of view that is being imaged, makes it possible to reconstruct the image. For the ultrasound waves used in medical imaging, which have wavelength on the order of a millimeter, an error of a millimeter in the positions of the transducers can ruin the quality of the reconstructed image. But it may not be practical to have the transducers rigidly fixed in position to within better than a millimeter, over the entire image acquisition time, when imaging a large and soft part of the body, such as the trunk. To overcome this problem, WO 2022/123464 to Axelrod and Chaviv, cited above, describes a wearable sensing garment for sensing ultrasound waves from a body part. The sensing garment includes a mesh of fabric adapted to be worn on at least a body part of a wearer, and a plurality of ultrasound trasducers, emitters and receivers, disposed on the mesh fabric, wherein each of the ultrasound receivers is adapted to detect ultrasound waves returned from at least a body part of the wearer and further adapted to convey detected ultrasound waves to a controller. The ultrasound receivers are adapted to collect signals emitted by a plurality of ultrasound emitters randomly or deliberately embedded into the wearable garment. Because the relative positions of the ultrasound transducers may vary from subject to subject, and during image acquisition for a given subject, they perform “simultaneous localization and mapping,” using FWI to simultaneously reconstruct an image and determine the positions of the transducers with respect to the body parts. Including the transducer positions as additional variables that the FWI is solving for may make the inversion problem even larger, and more likely to be unstable, if the field of view has a large number of grid points.

If a separate procedure is used to find the positions of the ultrasound transducers, before the image is reconstructed, then reconstructing the image will not be made more computationally expensive and potentially unstable. If acquisition of the imaging data takes a long enough time so that the positions of the transducers change significantly, for example by more than 2% or 5% or 10% or 20% or 50% of an ultrasound wavelength, during the data acquisition, then the procedure optionally finds the transducer positions as a function of time. The inventor has found that machine learning can be used to find the positions of the transducers, even as a function of time, to sufficient accuracy, using as input the imaging data itself. In particular, this can be done with machine learning using a convolutional neural network, trained on full-wave calculations of simulated ultrasound signal data in a body model.

It should be noted that it makes sense to find the positions of the transducers as a function of time, in order to accurately find ray paths and travel times for ultrasound waves when using prestack depth Kirchhoff migration to calculate angle gathers as a function of image reconstruction coordinates in the body, even if the motion of the body with respect to the image reconstruction coordinate system is ignored. That may be true because even an error of a fraction of a millimeter in the positions of the transducers can lead to major errors in the reconstructed image, but ignoring motion of the body over a few millimeters will likely only lead to blurring of the image over a few millimeters, which is not a major degradation of image quality if the resolution is only intended to be about 2 mm. But it may be true that, if the surface of the body moves by a few millimeters, then taking into account the change in positions of the transducers but not taking into account the motion of the body with respect to the image reconstruction coordinates may lead to major errors in calculated reflection angles for points within a few millimeters of the surface of the body, for nearby source transducers and receiver transducers. Optionally the reconstructed image is not found for locations close to the surface of the body, for example within a few millimeters of the surface of the body, or within 1 or 2 or 3 or 5 times the range of radial motion determined for the transducers.

In the methods and systems described herein, there are a plurality of source transducers and a plurality of receiver transducers mounted on the outside of a field of view in the body, and signals transmitted from the source transducers to the receiver transducers are recorded, and used to reconstruct an image of the field of view. Each recorded signal is transmitted by one of the source transducers to one of the receiver transducers. Optionally, only some of the possible combinations of source transducer and receiver transducer are used to generate a recorded signal. Optionally, there is a subset of source transducers that are used to transmit signals at all, and only a subset of receiver transducers are used to receive signals at all, but every source transducer in the subset of source transducers is used to transmit signals, that are recorded, to every receiver transducer in the subset of receiver transducers. For example, every other source transducer is used, and/or every other receiver transducer is used, to generate a lower resolution image than if all the transducers were used. Optionally, less than 20% of the source transducers, or between 20% and 50% of the source transducers, or between 50% and 80% of the source transducers, or more than 80% of the source transducers, are used to generate recorded signals. Optionally, less than 20% of the receiver transducers, or between 20% and 50% of the receiver transducers, or between 50% and 80% of the receiver transducers, or more than 80% of the receiver transducers, are used to generate recorded signals. Optionally, there is not a subset of the source transducers and a subset of the receiver transducers, such that every source transducer in the subset transmits a recorded signal to every receiver transducer in the subset, but there is a subset of source-receiver combinations that is used to generate recorded signals, for example less than 10% of the possible combinations, or between 10% and 20%, or between 20% and 50%, or between 50% and 80%, or between 80% and 90%, or more than 90% of the possible combinations are used to generate recorded signals. Optionally, a recorded signal is transmitted from every source transducer to every receiver transducer. This situation may allow the best quality images to be reconstructed. However, with the very large number of source and receiver transducers that may be used, it may be common for a few of them not to be functioning properly, and this may have relatively little effect on image quality. In the rest of the description of methods and systems herein, it will be assumed, for simplicity, that every source transducer transmits a recorded signal to every receiver transducer.

In many of the embodiments described herein, pre-stack depth Kirchhoff migration is used to generate an image from the recorded signals, using gathers, such as angle gathers, that are constructed for some points in the field of view. For example, there is a coordinate grid in the field of view, and gathers are constructed for some subset of the grid points in the coordinate grid, for example less than 10% of the grid points, or between 10% and 20% of the grid points, or between 20% and 50% of the grid points, or between 50% and 80% of the grid points, or between 80% and 90% of the grid points, or more than 90% of the grid points. Optionally, 100% of the grid points are used to construct the gathers. In the rest of the description of methods and systems herein, it will often be assumed, for simplicity, that all of the grid points in the field of view are used to construct the gathers. But in some embodiments of the invention, especially if the coordinate grid extends outside the body, gathers are not constructed for all grid points. An aspect of some embodiments of the invention concerns a system for generating ultrasound medical images of a field of view of a human or animal body, in which there are ultrasound transducers mounted in an array on the outer surface of the body on opposing sides of the field of view, and pre-stack depth Kirchhoff migration is used to generate the image from the data of ultrasound signals transmitted by source transducers and received and recorded by receiving transducers in the array. In Kirchhoff migration, information about each point in the field of view, for example each grid point in a coordinate system, is obtained by gathering the recorded signals received by each receiver transducer from each source transducer, evaluating each signal at a time delay after its transmission time by the source transducer, defined in some way, for example defined as the time of transmission of the beginning of the signal, that corresponds to the travel time for a wave emitted by that source transducer to reflect from that point and to travel to that receiver transducer. In the field of geological acoustic imaging, a recorded signal from a particular source and recorded by a particular receiver is often referred to as a “trace.” Although each of these signals also includes contributions from ultrasound waves that reflected from many other points in the field of view in going from that source transducer to that receiver transducer in that travel time, those contributions tend to cancel out by destructive interference when the signals for a given point are added up, and the contributions from ultrasound reflecting from that point tend to reinforce each other by constructive interference when the signals for that point are added up. In Kirchhoff migration, the travel time, and optionally the path of propagation from the source transducer to the point and from the point to the receiver transducer, may be found for example by ray tracing, or by directly solving the eikonal equation. In pre-stack migration, including in pre-stack Kirchhoff migration, all these signals are not simply scaled and added together for a given point, but they are sorted out by a parameter, such as reflection angle of the ultrasound waves at the point, or offset between the source transducer and receiver transducer, to form gathers, such as angle gathers or offset gathers. These gathers contain information that is optionally used to correct a velocity model of the field of view, giving the sound speed as a function of position in the field of view, starting with an initial velocity model. Since the velocity model is used to find the paths and travel times of the ultrasound waves when finding the gathers, which in turn are used to correct the velocity model, the correction of the velocity model can be done iteratively, until a convergence criterion is satisfied. Correcting the velocity model can be done even for systems and methods that have different features than the main embodiments described here.

The data in the gathers is also used to generate an image of the field of view, for example a reflectivity image, and/or an image of one or more other pre-stack attributes of the field of view, which in turn can be used to generate an image of a post-stack attribute, such as Poisson’s ratio or S-wave velocity. In depth pre-stack migration, such as depth pre-stack Kirchhoff migration, a velocity model is used to calculate ray paths for the ultrasound, which makes it possible to reconstruct the image accurately, as opposed to time migration, in which only the average velocity down to each depth is considered, and the ultrasound is assumed to propagate in straight lines, leading to easier computation but less accurate images.

Optionally the gathers are angle gathers, and include reflection angles greater than 50 degrees, or greater than 60 degrees, or greater than 70 degrees, or greater than 80 degrees. Optionally, at least some of the source transducers and receiver transducers whose signal data contributes to the high reflection angles for those angle gathers have a direct transmission path from the source to the receiver through the field of view. Optionally, the angle gathers exclude a range of reflection angles close to 90 degrees, for example within 10 degrees or within 5 degrees or within 2 degrees of 90 degrees, because it may be difficult to distinguish a signal that reflects from a point at an angle close to 90 degrees, from a signal that is directly transmitted from a source transducer to a receiver transducer without reflecting. Optionally, a reflection angle A that is greater than 90 degrees for a given point in the field of view is not included in that angle gather, but is treated as a reflection angle of 180° - A for an angle gather on the other side of the field of view, extending to the same point from the opposite direction.

Having angle gathers with a greater range of reflection angles has the potential advantage that there will be more signal data contributing to the image for a given point in the field of view, which may make images of the field of view more accurate and/or less noisy than if a smaller range of reflection angles were used. This may be especially true for images of attributes that depend on how the reflectivity depends on reflection angle, such as Poisson’s ratio. In addition, in ultrasound medical imaging it is possible to measure signals transmitted directly from a source to a receiver through the field of view without reflecting, or with only slight scattering, which may provide useful information about the field of view, for example information that could be used for ultrasound tomography.

It should be noted that it is generally possible to use a greater range of reflection angles for ultrasound medical imaging using angle gathers than for seismic imaging using angle gathers, because in ultrasound medical imaging source and receiver transducers may be located on multiple sides of the field of view, while in seismic imaging the acoustic sources and receivers are generally all located nearly in the same plane, on or slightly below ground level. To obtain data on large reflection angles for seismic imaging, for example for reflection angles greater than 60 or 70 or 80 degrees, for a point deep below ground level, would require using sources and receivers that are far away from the point horizontally in opposite directions, which would result in a weak signal, and would require installing geophones and receivers over a wide area, and this is generally not done. In addition, in seismic imaging, signals with large reflection angles often include a large amount of interference from “Head Waves,” and are not used for that reason. Head Waves are described, for example, by Vlastislav Cervey and Ravi Ravindra, Theory of Seismic Head Waves, University of Toronto Press, 1971, which may be accessed at <www(dot)jstor(dot)org(slash)stable(slash)10(dot)3138(slash)j(dot)cttlgxxrOg>. In seismic imaging using angle gathers, the maximum reflection angle used is typically no more than 40 degrees. But with ultrasound medical imaging using angle gathers, with transducers located all the way around the body part being imaged, or located a large fraction of the way around the body part, it is possible to obtain angle gather data for large reflection angles, using source and receiver transducers that are not too far away. This makes it possible to obtain data about the field of view in ultrasound medical imaging that would not be available, or would be much lower quality, in seismic imaging.

Optionally, the field of view is in a part of the body such as the trunk, or the limbs, and the transducers are mounted at different positions on opposing sides of the surface of the body part, for example by mounting the transducers on a flexible, stretchable frame, which is then wrapped around and fitted to the body part, which can have a different size and shape for different people. In these circumstances, the relative positions of the different transducers are not known in advance. In some embodiments of the invention the positions of the transducers can even vary over time, for example due to the subject breathing, while ultrasound signal data is being acquired for imaging the field of view.

Although pre- stack depth Kirchhoff migration, and other types of migration, are known for producing seismic images of underground geological formations, in these seismic imaging systems the transducers are fixed in place at different locations on the ground, and their exact locations can be measured, and do not change over time. Also, in these seismic imaging systems, the transducers are not located on opposing sides of the field of view, but are all located on the ground or buried near the surface, above all or almost all of the field of view, which is entirely underground. By contrast, in at least some embodiments of the medical ultrasound imaging system of the invention, transducers are located on all sides of the field of view, surrounding it. For points in most of the field of view, there is no plane passing through the point such that all of the transducers are located on one side of the plane. Whether the plane is oriented axially, coronally, or sagittally, or obliquely, in these embodiments there are always some transducers located on each side of it. Optionally, for definiteness, the location of a transducer is considered to be the location of the center of its face that is in acoustic contact with the body. In some embodiments of the invention, this is true separately for the source transducers and the receiver transducers — there are at least some source transducers located on each side of the plane, and at least some receiver transducers located on each side of the plane. As used herein, “location” of a transducer means the same thing as “position” of a transducer.

In some embodiments of the invention, not only are there transducers located on both sides of every plane that passes through a point in the field of view, for all points in most of the volume of the field of view, but this is also true for each subset of receiver transducers that records signal data from a given source transducer, for at least one source transducer, or for a plurality of source transducers, or for most source transducers, or for all source transducers. This is in contrast to some prior art ultrasound imaging systems where the field of view is divided into slices, and each slice is imaged using only ultrasound signals transmitted by source transducers located in or close to that slice to receiver transducers located in or close to that slice. For those systems, this condition is not satisfied, because, for points in most of the field of view, if we take a plane passing through that point that is parallel to the slice in which that point is located, and we consider any one of the great majority of source transducers that are located in different slices not very close to the slice of that point, then all of the receiver transducers that record signals from that source transducer will also be located in or near the slice where that source transducer is located, which will all be on the same side of the plane.

Conventional medical ultrasound imaging systems that use hand-held ultrasound transducer arrays, such as B-mode systems,, which do not use migration and other sophisticated techniques for reconstructing images, also have transducers that are located on only one side of the field of view, in a narrow area on the surface of the body. In these systems the transducers are arranged in a rigid array, with the relative positions of the different transducers fixed. The field of view is limited to a narrow region extending into the body from the transducer array, it is two-dimensional, and it can be translated or swiveled around by manually manipulating the instrument holding the transducer array. The person performing the exam views an image of the field of view in real time on a computer monitor, and decides when to save a screen shot of the monitor, using their medical expertise. This is in contrast to some, but not all, exemplary embodiments of the invention, where the field of view is three-dimensional and fixed, covering a large part of the volume of a body part that is being imaged, and the image is viewed not in real time, but later when it has been reconstructed. In those embodiments, acquisition of the imaging data is generally done automatically, optionally using a wearable transducer array, and optionally does not require medical expertise in real time. In some prior art medical ultrasound imaging systems that use migration, such as Theis and Bonomi, and Renaud et al, cited above, the ultrasound source and receiver transducers are all part of a conventional off-the-shelf rigid array of medical ultrasound transducers, located on only one side of the field of view. In such systems, the relative positions of all the transducers are fixed and known. That is also true of the ultrasound imaging system, used for non-destructive testing of concrete elements, described by Muller et al, cited above. In such prior art, transducers for ultrasound medical imaging would not be located on different sides of a non-rigid field of view, for example mounted on a soft body part, because their precise relative positions would not be known, so an image could not be reconstructed. In other prior art medical ultrasound imaging systems, the field of view is surrounded on multiple sides by ultrasound transducers, but the transducers are still held rigidly in place at fixed and known relative positions, and the image is reconstructed by methods other than pre-stack depth Kirchhoff migration, such as FWI or poststack migration and/or time migration. That is true of the published patent applications of Calderon Agudo, and the paper by Guasch et al, cited above, which describe systems for ultrasound imaging of the brain and skull using FWI to reconstruct the image. The transducers are mounted in a helmet that surrounds the head of the subject. Robins et al, cited above, describe such a system for imaging the breast, also using FWI, where the breast is immersed in water held in a rigid hemispherical container, with transducers mounted at fixed locations on the inner surface of the container. Ruiter et al, cited above, describes breast imaging in a similar configuration, with the image reconstructed by various methods of post-stack migration. Schmidt et al, cited above, also describes a system for imaging a breast immersed in a container of water, with a moving ring-shaped array of source and receiver transducers that only images one slice of the field of view at a time, using post-stack time migration to reconstruct the image. EP 1578274 to Duric et al, cited above, describes a system for imaging a breast compressed between two acoustic paddles, with moving 1-D transducer arrays that only image one slice of the field of view at a time, using post-stack migration to reconstruct the image.

There are also ultrasound medical imaging systems, with a field of view surrounded by transducers on multiple sides, in which the transducers are not rigidly held in place, but these systems also do not use pre-stack depth Kirchhoff migration to reconstruct the image. Axelrod and Chaviv, cited above, for example, describe ultrasound sensors disposed on the mesh fabric of a sensing garment worn on a body part, and use FWI for reconstructing the image as well as for determining the positions of the sensors with respect to the body part. Axelrod and Dafni, cited above, use a similar sensing garment, and describe using either FWI or various post-stack migration methods for reconstructing an image. They describe using a morphing function for dealing with small motions of the body.

An aspect of some embodiments of the invention concerns a medical ultrasound system for generating images of a field of view of a part of the body of a subject, using a flexible stretchable frame holding an array of ultrasound transducers, both source transducers and receiver transducers. Ultrasound signals are transmitted by the source transducers and received by the receiver transducers and recorded, and the recorded signals are used to reconstruct an image of the field of view. The frame wraps around the body part, adapting to the shape and size of the body part, which can vary among different people, and holding the transducers in good acoustic contact with the surface of the body part. For example, the frame is flexible and stretchable enough to adapt to different body parts, or to corresponding body parts of different subjects, that differ in one or more dimensions by at least 10%, or by at least 20%, or by at least 30%, or by at least 50%, or by at least a factor of 2, or by more than a factor of 2.

Because the transducers are not rigidly fixed in place relative to each other and to the body part, information is optionally obtained about positions of the transducers on the body of the subject, at the times that ultrasound signals are transmitted from source transducers to receiver transducers. Information on the positions of the transducers may be needed to accurately reconstruct an image from the recorded signals.

Precise information on the positions of the transducers is optionally obtained from the recorded ultrasound signals, after first obtaining initial estimates of the positions, for example by modelling the frame and the body, and optionally from sensors attached to the frame that measure its stretching and bending. Optionally, the precise information on the positions of the transducers is obtained from the recorded ultrasound signals before the recorded signals are used to reconstruct an image. This is in contrast to Axelrod and Chaviv, cited above, where FWI is used simultaneously to determine the positions of the transducers and to reconstruct an image. Optionally, the standard deviation of the error in the positions of the transducers is well below a wavelength of the ultrasound in the field of view, for example less than 1 mm or less than 0.5 mm or less than 0.2 mm or less than 0.1 mm. For example, this information is obtained using a computer programmed to run a program that is configured to improve its performance by training, such as a machine learning program, for example a neural net, or a program with any other configuration known to be useful for machine learning. The program is configured, for example as a result of training by a machine learning procedure in the past, to determine the positions of the transducers from the recorded ultrasound signals, for example to within an error with standard deviation less than 1 mm, or less than 0.5 mm, or less than 0.2 mm, or less than 0.1 mm. Optionally, the program is configured to improve its performance, by training by machine learning. For example, the program has a format that allows its performance to be improved by optimizing the values of a large number of parameters, such as weights in a neural network, by providing it with training examples and evaluating its performance for each training example, and calculating how to change the weights to further improve its performance Even if the weights have been completely optimized and no further improvement in performance is possible, the program may still be considered to be configured to improve its performance by machine learning, if the program has a format, such as a neural network, that would allow it to improve its performance if it had a different set of weights, not optimized.

Using a program trained by machine learning, such as a neural net, to determine the positions of the transducers may be especially suitable for an ultrasound imaging system that uses pre-stack depth Kirchhoff migration to reconstruct an image. That is because, in pre-stack depth Kirchhoff migration, in contrast to FWI, determining the positions of the transducers is done as a separate step, before reconstructing the image, and using a program trained by machine learning has been found by the inventors to be a good way to determine the positions of the transducers as a separate step.

Although the position information pertains to the times that ultrasound signals are transmitted from source to receiver transducers, obtaining the position information need not be done in real time, but may be done long after the ultrasound signal data is transmitted and received, using the recorded signal data, for example sometime before the image is reconstructed using the recorded signal data. Optionally, obtaining information about the positions of the transducers on the body of the subject comprises making estimates of the positions of the transducers on the body of the subject, for example using mechanical models of the subject’s body and of the frame of transducers, and information about the coordinates of the transducers on the frame, and optionally also using data about the mechanical state of the frame from strain sensors and/or angle sensors on the frame. Optionally, those estimates of positions of the transducers are used, together with the recorded signal data, as input to the program trained by machine learning, which then finds more precise positions of the transducers.

Here, “positions of the transducers” and “relative positions of the transducers” means the spatial positions of the transducers in place on the body, referring to the positions of the transducers relative to the body, or relative to the coordinate system used to reconstruct an image. As used herein, “spatial positions of the transducers” and “relative spatial positions of the transducers” have the same meanings, respectively, as “positions of the transducers” and “relative positions of the transducers.” These terms do not refer to the positions of the transducers in the coordinate system of the array, which remain unchanged even when the frame of the array is flexed or stretched on the body of the subject.

Optionally, in order to reconstruct accurate images, the positions of the transducers are determined within 2% or 5% or 10% or 20% or 50% of an ultrasound wavelength, which would be 1.54 mm for 1 MHz ultrasound at a typical sound speed of 1540 m/s in soft tissue. The inventor has found that is possible to achieve this accuracy with a program configured to be trained by machine learning, such as a convolutional neural network (CNN). The positions of the transducers need not be determined in real time, but may be determined later, before reconstructing the image from the signal data. The program optionally finds the positions of each pair of source and receiver transducers for which a signal transmitted from the source is recorded by the receiver, at the time when that signal is transmitted, even if the positions of the transducers change significantly during the data acquisition period. Optionally each source transducer only transmits a signal once, and it is recorded simultaneously by many receiver transducers. The acoustic travel time from the source to the receiver is typically less than a millisecond, too short for any significant change in positions of the transducers to occur, and the duration of each signal is typically even shorter. In referring herein to determining the position or relative position of a transducer at a particular time, the time when a signal is transmitted means essentially the same thing as the time when the signal is received. But the time needed to record signals from all the source transducers, one after the other, may be seconds or minutes, more than enough time for the positions of the transducers to change significantly.

Optionally, the transducers, mounted on the frame, also extend all the way around the body part, which has the potential advantage that it might be possible to obtain a better image that way. Here, “all the way around” means there is at least one path, encircling the body part, on which there is no gap between transducers that is greater than 20% of the path length, or no gap that is greater than 10% of the path length, or no gap that is greater than 5% of the path length, or no gap that is greater than 2% of the path length. Alternatively, the transducers only extend part of the way around the body part, covering only part of the frame, which has the potential advantage that the array may be less expensive, and the image reconstruction may require less computation, if there are fewer transducers. For example, the transducers extend at least 80% of the way around the body part, or at least 50% of the way around, or at least 30% of the way around, or at least 20% of the way around, with no gap greater than 10%, or no gap greater than 5%, or no gap greater than 2% of the distance around the body part, in the part of the frame that the transducers extend over.

Optionally, initial estimates for the positions of the transducers on the body of the subject are based on the position of each transducer in the coordinates of the flexible stretchable frame, and on a model of the body or body part of the subject that the frame is wrapped around, that estimates where each transducer will be located on the body of the subject. The model optionally includes body measurements of the body or body part of the subject, such as chest, waist and hip measurements if the frame is wrapped around the trunk of the subject, and/or one or more of the height, age and gender of the subject. Optionally, the frame comprises sensors, such as strain and angle sensors, that measure components of a two-dimensional strain tensor of the frame at different locations on the frame, which can be used to calculate more precise initial estimates for the positions of the transducers on the body of the subject.

In some embodiments of the invention, the program configured to be trained by machine learning is a neural network, such as a convolutional neural network (CNN).

Having an array of transducers that extends all the way or much of the way around a body part, for example the trunk, may make it possible to image a field of view that extends over all or much of the body part. For example, the array of transducers is wrapped around the trunk, and the field of view includes at least 50% of an axial cross-section of the trunk, or between 20% and 30%, or between 30% and 50%, or between 50% and 75%, or between 75% and 90%, or at least 90% of an axial cross-section of the trunk. Optionally, the field of view includes any of these percentages of the axial cross-section of the trunk if the axial cross-section is defined as only the area that includes internal organs and muscle, excluding skin and subcutaneous fat. Optionally, the array also extends some distance along the trunk axially, and the field of view is 3-D, extending axially along the trunk. For example, the field of view optionally includes all or most of one or more major organs, for example more than 50% of the volume of the heart and/or one or both lungs, and/or more than 50% of the volume of the stomach and/or the liver, and/or more than 50% of the volume of the small intestine and/or the colon, or between 20% and 30%, or between 30% and 50%, or between 50% and 75%, or between 75% and 90%, or at least 90% of one or more of these volumes. In this respect, this kind of ultrasound system can produce images that are more similar, in field of view, to the images produced by a CT or MRI system, than to the images produced by a conventional ultrasound medical system that uses a hand-held transducer array, such as a B-mode system. Because the images are reconstructed by migration, which reconstructs local information from waves that are received over a wide area, the quality of the images may also be closer to that of CT or MRI images, than to such conventional ultrasound images in which only the received wave data is displayed.

There are other potential advantages of using migration, especially pre-stack depth Kirchhoff migration, in medical ultrasound imaging with transducers surrounding the field of view, that do not occur with seismic imaging, where the transducer arrays extend only over a bounded area on the ground. As noted above, with transducers completely or largely surrounding the field of view, it is possible to obtain signal data for larger reflection angles, even close to 90 degrees or greater than 90 degrees, than in seismic imaging. In seismic imaging the transducers are all nearly in one plane corresponding to ground level, and reflection angles are typically limited to 40 degrees, and data for much greater reflection angles would require sources and receivers that are very distant from each other and would produce weak signals. Using data for a wider range of reflection angle may result in more accurate and/or less noisy images, especially for images of attributes that depend on reflection angle, such as density, or Poisson’s ratio. It is also possible, with transducers surrounding the field of view, to obtain data on direct transmission from sources to receivers without reflection. Such data could potentially be useful, for example for ultrasound tomography.

Furthermore, with seismic imaging, there are edge effects in a broad area near the edge of the transducer array, and that part of the field of view cannot be imaged accurately. In the case of a transducer array for medical ultrasound that extends all the way around the body part azimuthally, there are no edge effects in the azimuthal direction, although there may be edge effects in the axial direction, and more of the field of view can be accurately imaged. Another advantage to using prestack depth Kirchhoff migration with an array that goes all of the way or even part of the way around a curved surface of a body part, is that the gathers, such as angle gathers or offset gathers, may extend below the surface in a direction at least approximately normal to the surface, and, as a result of the curvature of the surface, the gathers may get closer together as they extend further beneath the surface. When the velocity model is updated, using the data of the gathers, information about the interval velocity at neighboring gathers at the same depth may be used to get a more accurate estimate of the interval velocity, for example by smoothing the interval velocity over the azimuthal coordinate. Furthermore, at a value of r far enough below the surface, where the gathers get close enough together in the azimuthal direction, some of the grid points can be removed from the coordinate grid for values of r higher than that value of r, for example grid points can be removed for some values θ, and the corresponding gathers can end at that value of r. Using fewer grid points can make the computation more efficient. These advantages are not available for seismic imaging when the ground is flat and gathers remain the same distance apart as they extend beneath the surface.

An alternative method to migration, often used in seismic imaging, is full waveform inversion (FWI), in which the local sound speed as a function of position in the field of view is calculated from the recorded seismic signals by inverting the solution to the full wave equation. Although FWI is computationally expensive, and, unlike migration, it requires a fairly accurate initial velocity model of the field of view, it has the advantage over migration in seismic imaging that it is not necessary to manually adjust different parameters depending on characteristics of the geological formation being imaged. Calderon- Agudo, cited above, describes a system for medical ultrasound imaging of the head, in which transducers are rigidly mounted in fixed relative positions on a helmet that goes over the head, and the image is reconstructed from the ultrasound signal data using FWI. Because the helmet is rigid, and the subject’s skull is rigid, the transducers have known relative positions that do not change. In addition, it is possible to determine a fairly accurate initial velocity model of the field of view, because the speed of sound is much greater in the skull than in the brain, and the distribution of bone and soft brain tissue inside the head is approximately known in advance. But such a medical ultrasound system, with transducers mounted around the field of view, could not be easily adapted for imaging the trunk or other body parts that do not have a rigid outer surface, and where the speed of sound does not vary much over the field of view. In addition, FWI would be especially computationally expensive, and might not be practical for imaging a large field of view such as a large part of the trunk at high resolution, also because of the possibility that the computation might be unstable. And the need to do the calculations over again for every new configuration of the transducers, as the soft outer surface of the trunk changes shape, would make FWI even more computationally expensive.

An aspect of some embodiments of the invention concerns medical ultrasound imaging of a field of view of the body, using automated workflow. In seismic imaging of geological formations, it is generally necessary to manually select the various parameters used for reconstructing the image with pre-stack depth Kirchhoff migration, because geological formations can greatly differ from each other. Extensive manual intervention is especially needed for the process of correcting the velocity model using angle gathers or offset gathers of the signal data. Examples of this need for manual decision making are given in the paper by Bidikhova and Hall, cited above. The inventor has found that it is possible to avoid this problem in ultrasound medical imaging, when using pre-stack depth Kirchhoff migration, because human bodies are much more similar to each other, in how sound speed varies spatially, than geological formations, because human bodies all have the same organs, of nearly the same shapes at nearly the same positions in the body. Once a data acquisition device has been designed, and a data collection protocol and other software parameters have been carefully selected that gives good results for imaging a given body part of one subject, it will generally work satisfactorily for imaging the same body part of other subjects. This allows fully automated medical ultrasound imaging using pre-stack depth Kirchhoff migration. The hardware parameters for the data acquisition device include the number of transducers, their locations in the array, the shape, topology and elastic properties of the array, and the peak frequency and bandwidth of the transducers. The parameters that are specified, including software parameters, for reconstructing a medical ultrasound image from data of the received ultrasound signals using pre-stack depth Kirchhoff migration, include, for example, an initial velocity model of the field of view; the resolution and shape of the coordinate grid for the field of view; the sampling rate of the signal data; how many and which transducers are used as source transducers, and the signal envelope they use; parameters for how the recorded data of received ultrasound signals is used to create gathers, such as angle gathers or offset gathers, for example the range and binning of reflection angles used for angle gathers; the parameters used for correcting the velocity model using the data in the gathers, for example the number of iterations and/or other convergence criteria for correcting the velocity model; and the parameters used for finding pre-stack attributes from the gathers, and the parameters used for finding post-stack attributes from the pre-stack attributes.

For purposes of better understanding some embodiments of the present invention, as illustrated in FIGS. 2 to 13 of the drawings, reference is first made to the construction and operation of a prior art seismic imaging system, as illustrated in FIG. 1. Figure 1 shows a prior art seismic imaging system 100, with an array 102 of acoustic sources, and geophones that receive sound waves generated by the acoustic sources and reflected by underground features. System 100 uses the signals generated and received by array 102 to reconstruct an image of an underground field of view 104. Individual sources and geophones 106 in the array are arranged on the surface of the ground above field of view 104, typically in a Cartesian grid, and are fixed in place on the ground, in good acoustic contact with the ground, or buried a short distance under the ground. Several publications reviewing methods of seismic acoustic imaging are cited above.

A number of methods are available for reconstructing seismic images from the signals, including full waveform inversion (FWI), where the full 3-D wave equation for sound waves in anisotropic media, for a given sound speed and anisotropy as a function of position in the field of view, is inverted to produce a 3-D map of the sound speed and anisotropy. FWI requires that an approximate map of the sound speed in the field of view is known in advance. It is very expensive computationally. Because even a single seismic image of an underground field of view may contain valuable information about oil and other minerals, and because buying and placing the sources and geophones may also be very expensive, it may be cost effective to use FWI to reconstruct the image in spite of the cost of computation.

Various types of migration, including Kirchhoff migration, may also be used to reconstruct the image, and these methods are typically less computationally expensive than FWI, and do not necessarily require having an approximate map of the sound speed in advance, since even a uniform sound speed within the field of view can be used as an initial velocity model. However, migration methods typically have a lot of free parameters that have to be chosen manually depending on geological characteristics of the field of view, and these methods require highly trained specialists to perform them.

Before explaining at least one embodiment of the invention in detail, it is to be understood that the invention is not necessarily limited in its application to the details of construction and the arrangement of the components and/or methods set forth in the following description and/or illustrated in the drawings and/or the Examples. The invention is capable of other embodiments or of being practiced or carried out in various ways.

Medical ultrasound imaging system

Referring now to the drawings, FIG. 2 illustrates an ultrasound imaging system 200 comprising an array 202 of ultrasound transducers 208 that are in good acoustic contact with an outer surface of a body part of a subject. The transducers comprise both source transducers that transmit ultrasound waves into the body, and receiver transducers that receive ultrasound signals that have reflected from structures inside the body, as well as ultrasound waves that are transmitted directly through the body from source transducers. Optionally, the source transducers and the receiver transducers have identical construction, and optionally some or all transducers are used some of the time as source transducers and some of the time as receiver transducers, during an imaging procedure.

Optionally, array 202 is flexible and stretchable, and can be adjusted to fit the outer surface of the body of many different subjects, whose bodies or body parts have a range of different sizes and shapes, and to hold the transducers in good contact with the body surface. Optionally, array 202 is configured, for example by size, shape, structural arrangement, elastic properties, and/or fastener configuration, to wrap completely or part way around the body part being imaged, for example the trunk, and system 200 images a field of view 204 inside the body part, with transducers located on more than one side of the field of view. Array 202 is optionally constructed of a flexible and stretchable frame 206, for example with a grid-like structure, with individual transducers 208 mounted on the frame, for examples at the intersections of the grid. The grid optionally is rectangular, Cartesian, hexagonal, triangular, or has any other regular or irregular geometry. The grid lines may be numbered, or otherwise labeled, as a way of specifying coordinates where each transducer is located on the frame. For example, in the case of a Cartesian grid, the grid lines going in one direction, for example vertically in FIG. 2, may be referred to the warp, and the grid lines going the other way, for example horizontally in FIG. 2, may be referred to as the woof, and the position of any transducer on the frame may be defined by its warp and woof coordinates. A different coordinate system may also be used, for a Cartesian grid or any other type of grid, and the coordinate system need not correspond to numbered grid lines, but could be, for example, a transformation of numbered grid lines. The coordinates that specify the location of a transducer on the frame are referred to herein as its frame coordinates, or array coordinates, to distinguish these coordinates from the spatial coordinates of the transducer when the frame is in place on the body of the subject. Alternatively, frame 206 of array 202 is constructed from a continuous sheet of elastic material or a woven elastic cloth, rather than having discrete grid lines as in FIG. 2, with transducers mounted at certain locations on the sheet. In any case a coordinate system may be defined on the sheet, or on any other construction of the frame of the array, in order to specify the locations of transducers in the array, and these coordinates are referred to herein as frame coordinates or array coordinates. Optionally, the array has a design specification that assigns a unique ID number to each transducer, and lists the coordinates of each transducer on the array, together with its ID number. The design specification of the array may also specify an elastic model of the array, which can be used to estimate how the transducers will be distributed on the surface of the subject’s body when the array is in place on the body.

The array, before it is wrapped around the trunk, or another approximately cylindrical body part such as a limb, optionally has a two-dimensional shape that can be stretched to fit snugly around the trunks of subjects with a wide range of body shapes, for example subjects whose body measurements, for example chest, waist and hip measurements in the case of an array that goes around the trunk, fall within certain relatively broad ranges. Optionally, the array is made in different sizes to fit subjects with body measurements that fall within different ranges.

Optionally, the array has a mechanical structure, for example an anisotropic elastic structure, that will cause it to conform closely to the body surface of the subject, including parts of the body surface with negative curvature or concavities, where a simple isotropic elastic sheet might tend to pull away from the body surface. For example, the array has a structure like the structure of a mechanical counterpressure (MCP) suit, which was the subject of research in the 1950s and 1960s by NASA, which hoped to develop such a skin-fitting mechanical pressure suit for astronauts as an alternative to the bulky inflated air pressure suits that NASA ended up adopting for astronauts. This research, as well as ongoing research on MCP suits at MIT, is described in the Wikipedia article on “Mechanical counterpressure suit, <en(dot)wikipedia(dot)org/wiki/Mechanical_counterpressure_suit>, downloaded on September 28, 2022. Optionally, a flexible stretchable frame for an ultrasound array does not exert pressure, or does not exert adequate pressure, on one or more problematic areas of the body surface, for example small concave or negative curvature areas of the body surface, and transducers are absent from those areas, or any transducers on those areas are effectively absent because they are not pressed into good acoustic contact with the body surface. The array may still produce a good image using data from transducers located in other areas, where the transducers are in good acoustic contact with the body surface. Alternatively, transducers in the problematic areas are held in good acoustic contact with the body surface by an adhesive layer that transmits ultrasound, for example an adhesive layer similar to the adhesive layer described in the MIT News article cited above.

One or more fasteners 210 optionally hold the ends of the frame together when it is wrapped around the trunk or body part. The fasteners may be, for example, hooks, or buttons, or Velcro®, or a zipper. Optionally, sensors 211, for example strain sensors and/or angle sensors, are incorporated into frame 206, to measure distances and/or angles between adjacent transducers when the array is in place, stretched over the surface of the body. Information from these sensors may be used by the image reconstruction algorithm to determine the positions of the transducers in real time, as will be explained below.

It should be noted that FIG. 2 is drawn in a highly schematic way and is not at all to scale. Typically, the size and spacing of the transducers will be comparable to the desired resolution of the image, and the wavelength of the ultrasound waves used. For example the transducers are spaced, in one or both directions along the array, at intervals of less than 1 mm, or between 1 and 2 mm, or between 2 and 5 mm, or between 5 and 1 cm, or between 1 and 2 cm, or more than 2 cm, and the diameter of each transducer is optionally a significant fraction of the spacing, for example the transducers are less than 10%, or between 10% and 20%, or between 20% and 30%, or between 30% and 50%, or between 50% and 70%, or between 70% and 90%, or more than 90%, of the spacing. Optionally, the transducers have dimensions that make them resonant at the ultrasound frequency that is used for the signals, for example less than 200 kHz, or between 20θ and 500 kHz, or between 500 kHz and 1 MHz, or between 1 and 2 MHz, or more than 2 MHz. Using a lower frequency has the potential advantage that the ultrasound waves are less damped in body tissue and penetrate deeper into the body, but using a higher frequency has the potential advantage that the image can have a higher resolution. Optionally the frame is adapted to wrap snugly around a trunk with a maximum circumference that is less than 75 cm, or between 75 and 100 cm, or between 100 and 125 cm, or between 125 and 150 cm, or between 15θ and 175 cm, or greater than 175 cm.

Optionally, the transducers extend over all or almost all of the frame, and are densely spaced, for example at a spacing in one of the ranges mentioned above, over the whole body part surface, or almost all of the body part surface, for example more than 90% or more than 95% or more than 98% of the body part surface, that the frame is wrapped around, when the frame is wrapped around the body part and fastened. This has the potential advantage that transducers will be located on all sides of the field of view, at least azimuthally around a body part, such as the trunk, that is approximately cylindrical. Optionally, the array, and the field of view, also extends some distance axially along the body part, for example by less than 20% of an average diameter of the body part where the array surrounds it, or between 20% and 50% of the average diameter, or between 50% and 100% of the average diameter, or between 1 and 2 times the average diameter, or more than 2 times the average diameter. Optionally, for all points in most of the volume of the field of view, there is no plane passing through the point, such that all of the transducers are on one side of the plane, and none of the transducers are on the other side. That would be true, for example, of array 202 and field of view 204 in FIG. 2. As discussed above, there are advantages to having transducers surrounding the field of view, for example the absence of edge effects, at least in the azimuthal direction.

Alternatively, the transducers only cover part of the body surface azimuthally when the frame is wrapped around the body part, for example less than 30% of the way around, or between 30% and 50% of the way around, or between 50% and 70% of the way around, or between 70% and 80% of the way around, or between 80% and 90% of the way around, or more than 90%. Even if the transducers only extend part of the way around the body part azimuthally, optionally the frame itself still goes all the way around the body part, so that its edges can be fastened together at one or more axial locations, and its tension holds the transducers snugly against the body surface, in good acoustic contact. Having the transducers go only part of the way around the body surface has the potential advantage that, with fewer transducers, less computation may be needed to reconstruct the image, and even if there are edge effects, the edge effects may only affect a part of the field of view that is not important for the medical purpose of the imaging procedure.

A controller 212, for example a computer, or dedicated circuitry, controls the signals transmitted by the source transducers, through a communications connection 214, such as a cable, and records the signals received by the receiver transducers, through a communications connection 216, such as a cable. The controller controls the shape and amplitude of the waveform transmitted by each source transducer, and controls the timing of transmission by each source transducer, as well as recording the timing of signals received by the receiver transducers, and optionally digitizing the signals for storage as data. Alternatively, some of those functions are performed by elements in the transducers themselves, or incorporated into the array, but in that case those elements may be considered to be part of the controller. The controller need not be a single unit, such as a single computer, but its functions are optionally distributed among a plurality of units. Optionally, during the imaging procedure, only one source transducer transmits a signal at a time. Optionally, for each transmission by a source transducer, the signals are recorded from a plurality of receiver transducers, optionally from all of the receiver transducers, which optionally comprise all of the transducers that are not acting as the source transducer. The signal from each receiver transducer is recorded by the controller, with the data of each signal associated with an identifier of the source transducer that transmitted the signal, and an identifier of the receiver transducer that received the signal. The controller can distinguish the signals received from each receiver transducer, for example because different receiver transducers communicate with the controller over different cables or communication channels, or because different receiver transducers operate at different times, during the transmission of the signal by one source transducer, so it can associate an identifier of the receiver transducer and an identifier of the source transducer with that data.

Controller 212, or optionally a separate computer or controller which will be considered to be part of controller 212, analyzes the stored signal data to reconstruct one or more images of the field of view, using an image reconstruction algorithm that will be described. With this algorithm, the images optionally include a main image, for example showing reflectivity of the body tissue in the field of view, and one or more additional images showing one or more attributes of the field of view, such as sound speed, or Poisson’s ratio, as a function of position in the field of view. The image reconstruction algorithm optionally depends on knowing the position of the source transducer and the receiver transducer for each piece of the recorded signal data.

A power supply 218 supplies power for controller 212, as well as supplying power to the source transducers for transmitting ultrasound waves, and power to the receiver transducers for producing an analog signal, and/or digital data, of the received signal, and sending it to the controller for recording. Optionally, the power to the transducers is delivered from controller 212 through connections 214 and/or 216. Alternatively, some or all of the transducers are powered directly from power supply 218. Optionally, power supply 218 is not a single physical unit, but represents a plurality of physically separate power supplies. A display 220 optionally displays the image or images reconstructed by controller 212, in real time, or later, when requested by medical personnel who are using the image or images. In contrast to conventional medical ultrasound systems using a hand-held ultrasound array, where a 2-D image is almost instantaneously visible on a display screen to the person doing the procedure, the images reconstructed by system 20θ are much better quality, and may be fully 3-D, but generally take much longer to process, and may only be available sometime after the procedure is performed. Figure 3 shows a flowchart 300 for a method of making a medical ultrasound image of a subject, according to an exemplary embodiment of the invention. At 302, a subject is provided, and the body part to be imaged is specified. The body part could be, for example, the trunk, with the field of view covering all or most of a cross-section of the trunk, extending axially over a substantial part of the volume of the trunk, or the body part could be another cylindrical body part, such as an arm or leg, with the field of view covering all or almost all of a cross-section of the arm or leg, and extending some distance along the arm or leg. This is in contrast to conventional ultrasound imaging systems using a hand-held ultrasound transducer array, where the field of view typically only extends a short distance in from the surface of the body.

At 304, an ultrasound array with a stretchable, flexible frame, such as array 202 in FIG. 2, is wrapped around the body part to be imaged, and fastened to hold it in place, and to hold the transducers of the array in good acoustic contact with the outer surface of the body part. A gel is optionally used to facilitate transmission of ultrasound waves between the source and receiver transducers and the body tissue, as is typically done when making ultrasound images. Alternatively, a rigid ultrasound array is used, and it is pressed against the surface of the body rather than wrapped around it and fastened. But using a flexible stretchable array, wrapped all the way around, has the potential advantages that it can conform well to the body shape of many different subjects, that it doesn’t have to be held in place during the imaging process, that the elastic tension of the array can exert a force holding the transducers against the body, and that having transducers going all the way around the body part, or a large fraction of the way around the body part, potentially can produce a better quality image covering a large and 3-D field of view.

At 306, parameters that will be used for the image reconstruction are provided. The parameters include an initial velocity model of the body part, migration parameters that are used for creating gathers of the signal data, such as angle gathers or offset gathers, that can be used for correcting the velocity model, and parameters that are used for correcting the velocity model, as well as parameters for generating images from the data of the gathers. Examples of migration parameters, parameters for correcting the velocity model, and parameters for generating images, include the resolution and shape of the coordinate grid, the sampling rate of the signal data, the range and binning of reflection angles used for angle gathers, and the number of iterations and/or other convergence criteria for correcting the velocity model.

Optionally, especially if the field of view is mostly soft tissue, the initial velocity model sets the sound velocity at all locations equal to the speed of sound in water, 1540 m/s at body temperature and at the salinity of blood, which is generally a very good approximation for soft tissue. Even if a small part of the field of view comprises bones, such as ribs, hips, and spine in the case of a trunk image, or the interior bones of the limbs in the case of a limb image, the process of correcting the velocity model will generally converge at a sufficiently accurate sound speed for most of the field of view so that an accurate reflectivity image can be generated using migration.

In contrast to the situation in seismic imaging, where migration parameters have to be chosen manually by experts depending on the type of geological formation being imaged, with medical ultrasound imaging, at least for a given body part, all human subjects are pretty similar, and a single set of parameters, worked out in advance, can be often be used for all subjects, and can give good results. Furthermore, because of the similarity in human bodies, the same ultrasound imaging device, using the same protocols, can be used for many different subjects. Optionally, some adjustments in parameters are made automatically depending on characteristics of the subject that are known in advance, for example age, weight, height and gender, and it is not necessary for an expert to manually choose migration parameters. The parameters also optionally depend on the type of image being made, for example which body part is being imaged, the required resolution, and the type of contrast agent used if any. In addition, the parameters are optionally adjusted in response to information obtained from the transducer array when it is wrapped around the body part, for example information from strain sensors in the frame, about how much and where the frame is stretched.

At 308, ultrasound signals are transmitted into the field of view of the body part, optionally from one source transducer at a time. This is a major embodiment of the invention that is presently being considered by the inventor, and it has the potential advantage, in the usual case where the width of a transducer is comparable to or smaller than the ultrasound wavelength, that the transducer transmits ultrasound signals to a large number of receiver transducers over a broad range of directions simultaneously, and these signals can all be recorded simultaneously. In some other embodiments of the invention also being considered by the inventor, a plurality of source transducers transmits signals simultaneously, for example a set of source transducers within a small distance of a central source transducer, for example within 2 mm, or within 5 mm, or within 10 mm of a central source transducer. This has the potential advantage that the different source transducers can be phased to produce a transmitted ultrasound beam with a limited angular spread, and that information can be taken into account by the migration algorithm to produce a more accurate image, as will be described below for FIG. 4. Optionally, in this case, more than one transmission is made from the plurality of source transducers at different times, with each transmission having a different phase relation between the transducers, directing the transmitted ultrasound beam in a different direction. Transmitting ultrasound from a plurality of source transducers together, with a well-defined phase relation between them, has the potential advantage of producing a higher signal-to-noise ratio for the received signal. Also, if the source transducers are all assigned to groups and each source transducer is only used together with the other source transducers in its group, this has the potential advantage that the total number of recorded signals may be lower than if each source transducer transmitted separately, so recording signal data may take less time and calculating images may take less time, but the total transmitted ultrasound power and lateral resolution may still be the same.

Optionally, the transmission of signals at 308 is gated to the breathing cycle and/or to the cardiac cycle of the subject, which has the potential advantage that it may reduce motion artifacts in the reconstructed image due to breathing or blood circulation. Optionally, in this case, different images are reconstructed for different phases of the breathing cycle and/or different phases of the cardiac cycle. Additionally or alternatively, signals are transmitted only during time intervals when the subject is holding their breath, and there may be several such intervals, spread out in time, so that there is enough total data acquisition time to acquire all the signal data needed to reconstruct an image.

At 310, for the transmission from each source transducer, the signal received by each of a plurality of receiver transducers, being used with that source transducer, is recorded as signal data, associated with identifiers for the source transducer and the receiver transducer. The signal data is recorded at time intervals 6t between samples, where 6t will be referred to as the sampling interval. For example, the same time interval 6t is used between all successive data samples, or 6t is the longest time interval between successive data samples. Optionally, the signal data is recorded with a sampling interval that is at least several times shorter than the central ultrasound wave period, for example at least 2 times or at least 5 times or at least 10 times or at least 20 times shorter than the central ultrasound wave period, or the wave period of the maximum significant ultrasound frequency component, in order to accurately preserve information about the phase of the signal, and to avoid aliasing effects associated with the sampling interval. Accurately preserving information about the phase of the signal is potentially important, because of the role that constructive and destructive interference plays when using the signal data for Kirchhoff migration, as will be described below in FIG. 4.

Optionally, all of the receiver transducers are used for each source transducer. Alternatively, only some receiver transducers are used for each source transducer, for example only receiver transducers that are located within a specified distance of the source transducer, for example within 10 cm or within 20 cm or within 30 cm. Using only some receiver transducers for each source transducer, or using fewer receiver transducers, has the potential advantage of reducing the amount of signal data that is stored and reducing the processing time for generating images. But using all receiver transducers, or using more receiver transducers, for each source transducer, has the potential advantage of producing a higher resolution image and/or a more accurate image. Also, using at least some receiver transducers that are further away from the source transducer has the potential advantage of allowing angle gather data that includes larger reflection angles, especially for locations that are deeper inside the body, which may produce more or better data about pre-stack attributes of the body tissue.

Optionally, the other source transducers, that are not being used to transmit the ultrasound, also serve as receiver transducers for the source transducer that is transmitting the ultrasound. Alternatively, some or all of the source transducers are never used as receiver transducers, and optionally source transducers that are never used as receiver transducers have a different structure than transducers that are used only as receiver transducers or that are used both as source and receiver transducers at different times, for example in order to optimize their performance as source transducers. Optionally, only a relatively small fraction of the transducers, for example chosen randomly, and/or distributed uniformly, are ever used as source transducers. For example, between 1 in 2 and 1 in 5, or between 1 in 5 and 1 in 10, or between 1 in 1θ and 1 in 100, or between 1 in 100 and 1 in 1000, or less than 1 in 1000, of transducers are used as source transducers. Alternatively, all of the transducers, or a substantial fraction of them, for example more than half of them, are used as source transducers, for example one after the other. Some potential advantages of using only a small fraction of the transducers as source transducers, especially if there are a large number of transducers, as there may be particularly if a large body part such as the trunk is being imaged at high resolution, is that it will take less time to obtain data from all the source transducers, it will require less memory to store all the data, and it will require less computation to reconstruct the image. Using only a small fraction of the transducers as source transducers may not degrade the image very much compared to using all or a large fraction of the transducers as source transducers, because the data from all the source transducers may be largely redundant, and the transducers that are used may provide more than enough data to generate good images. Optionally, an interpolation algorithm is used to reconstruct what the signals would be if they were transmitted by one or more transducers that are not in fact used as source transducers, or what the signals, transmitted by a given source transducer, would be if they were received by one or more transducers that are not in fact used as receiver transducers for that source transducer. For example, using such an interpolation algorithm may produce images that are about as good as would be produced using about twice as many source transducers without an interpolation algorithm, or using about twice as many receiver transducers without an interpolation algorithm. An example of such an interpolation algorithm is the so-called 5-D interpolation algorithm, described for example by Daniel Trad, “Five-dimensional interpolation: New directions and challenges,” CSEG Recorder, vol. 39, March 2014, downloaded from <csegrecorder(dot)com(slash)view(slash)five-dimensional-interpolation- new-directions-and-challenges> on October 10, 2022, and by Satinder Chopra and Kurt J. Marfurt, “Preconditioning seismic data with 5D interpolation for computing geometric attributes,” SEG Houston 2013 annual meeting, pages 1368-1373, downloaded on October 10, 2022 from <dx(dot)doi(dot)org(slash) 10.1190(slash)segam2013(slash)2013-0169.1 > .

At 312, a migration algorithm is used to reconstruct one or more images of the field of view, from the recorded signal data received by each receiver transducer from each source transducer. Details of how this is done will be given in FIG. 4. At 314 the one or more images are outputted.

Migration and angle gathers

Figure 4 shows a flowchart 400, for a method of reconstructing one or more medical ultrasound images from ultrasound signal data received by each receiver transducer from each source transducer, using migration. Migration is a method for using signals recorded from many different combinations of source and receiver transducers to image the locations in a field of view where the signals were reflected, which can be used to create a reflectivity image of the field of view. In migration methods, a reflectivity image is found directly from the signal data, for example by putting all the signal data into a formula. This is in contrast to full wave inversion, where a wave velocity or impedance field as a function of location in the field of view is found from all the recorded signal data, and from all transmissions by source transducers, by using an inversion algorithm, starting with a trial wave field and solving an optimization problem iteratively.

There are a number of different migration methods known from seismic imaging. Some of them are described, for example, in the books by Biondi, by Berkhout, and by Robein, cited above. The inventor has determined that pre-stack depth Kirchhoff migration, including correction of the velocity model, is the most suitable migration method to use for ultrasound medical imaging using the system of FIG. 2 and the method of FIG. 3, because it is the easiest to adapt to a very large number of transducers with irregular geometry.

In pre- stack depth Kirchhoff migration, for each combination of source and receiver transducers, and for each target location in the field of view, a sound velocity model for the field of view is used to find the travel time, and optionally the ray optics path, for ultrasound waves transmitted from the source transducer, reflecting from the target location, and reaching the receiver transducer. To find the reflectivity of a given target location, all the recorded signals from any source-receiver pair in the data, that could have reflected from that target location, evaluating each signal at the time it would take to travel from that source to that target location and then to that receiver, are added together.

In order to cover all target locations in the field of view, the signal data for each sourcereceiver pair is optionally recorded for a time equal to twice the time an ultrasound wave could take to travel across the longest dimension of the field of view, for example 1 milllisecond if the field of view includes most of the trunk. Alternatively, each signal is recorded for less than 0.1 milliseconds, or between 0.1 and 0.2 milliseconds, or between 0.2 and 0.5 milliseconds, or between 0.5 and 1 millisecond, or between 1 and 2 milliseconds, or more than 2 milliseconds. If the sampling time is 0.1 microseconds, equal to 10% of the wave period of a 1 MHz ultrasound wave, and the signal is recorded for 1 millisecond, then there will be 10,000 data samples recorded for each signal. Alternatively there are less than 1000, or between 100θ and 2000, or between 200θ and 5000, or between 500θ and 10,000, or between 10,00θ and 20,000, or between 20,00θ and 50,000, or between 50,00θ a nd 100,000, or more than 100,000 data samples recorded for each signal. Optionally, some signals are recorded for a longer time, and have more data samples than other signals. Optionally, in that case, the shorter data signals are filled in, for example with zeroes, to make all the data samples the same length. The time data for each recorded signal is defined relative to a specified time in the transmission of the signal from the source transducer, for example at the beginning of the transmission of the signal.

Contributions to the signals due to reflection from that target location will tend to interfere constructively, while contributions due to other ray paths will tend to interfere destructively and cancel out. The method as described in flowchart 400 uses pre-stack depth Kirchhoff migration, including correction of the velocity model, although other migration methods may also be used within the bounds of the invention. But many other migration methods were developed for use in seismic imaging with sources and receivers arranged in on a flat surface, and involve solving a wave equation in a Cartesian coordinate system. Using pre-stack depth Kirchhoff migration has the potential advantage over using other migration methods, that pre-stack depth Kirchhoff migration may be better adapted to medical ultrasound imaging using transducers arranged on a surface of the body with irregular curved geometry and possibly cylindrical topology, using an arbitrary curvilinear coordinate system. Another potential advantage of using pre-stack depth Kirchhoff migration is that it makes it easier to construct angle gathers, which may be useful for determining certain attributes of the body tissue in the field of view, such as Poisson’s ratio.

At 404, an angle gather is generated for each of a set of discrete locations, designated by a pair of coordinates (θ, z), on the surface of the body. The number of angle gathers used is great enough to produce an image of a specified resolution, consistent with the number of source and receiver transducers used, and the number of signals used. For example, there are less than 1θ angle gathers, or between 1θ and 2θ angle gathers, or between 2θ and 5θ angle gathers, or between 50 and 10θ angle gathers, or between 10θ and 20θ angle gathers, or between 20θ and 50θ angle gathers, or between 50θ and 100θ angle gathers, or between 100θ and 200θ angle gathers, or between 200θ and 500θ angle gathers, or between 500θ and 10,00θ angle gathers, or more than 10,00θ angle gathers. If fewer angle gathers are used, the image may have lower resolution, even if the number of transducers and signals used would be adequate to support a higher resolution. If more angle gathers are used, the computation needed to generate the image may be higher. Optionally, the number of gathers is comparable to the number of transducers, for example the number of receiver transducers. For example, the number of angle gathers per area of the array is between 0.5 and 1 times the number of transducers per area of the array, or between 1 and 2 times the number of transducers per area, or between 2 and 5 times the number of transducers per area, or between 5 and 20 times the number of transducers per area, or less than 0.5 times or more than 10 times the number of transducers per area. Optionally, the angle gathers are fairly uniformly distributed over the outer surface of the body over the field of view. Optionally the number of transducers per area of the outer surface of the body over the field of view varies from place to place, for example it is greater on a part of the field of view where a higher resolution image will be obtained, and in that case the number of angle gathers per area also optionally varies from place to place, for example at least approximately in proportion to the number of transducers per area, which may produce the highest resolution image that can be obtained in each part of the field of view with the signal data available for that part of the field of view.

The angle gather gives the migrated ultrasound reflectivity M of an image location (r, θ, z) in the field of view, as a function of reflection angle a of ultrasound waves from the image location. To find M(r, θ, z, a) for a given image point (r, θ, z) and reflection angle a, use is made of all of the pairs (S, R) of source transducers S and receiver transducers R that are positioned on the surface of the body so that an ultrasound wave transmitted from source S will be received by receiver R after reflecting from image point (r, θ, z) with reflection angle a, for example using ray optics to find the path of the ultrasound waves. It should be noted that, in the context of the ray optics picture being used here, the transmitted wave is referred to as “reflecting” from a single image point (r, θ, z), but physically the scattering of the wave, with its finite wavelength, is properly referred to as “diffracting” from image point (r, θ, z), and this method is referred to in the literature on Kirchhoff migration as “diffraction summation.” Reflection angle is defined herein as half the angle between the ray path of the transmitted ultrasound and the ray path of the reflected ultrasound, where they join at image point (r, θ, z). A velocity model that gives the sound speed as a function of position inside the body is used to calculate the ray optics paths. Also, for each source-receiver pair (S, R), the velocity model is used to calculate the total ray optics travel time TS,R(T, θ, z) from the source S, to the image point (r, θ, z), to the receiver R.

Figures 5 A and 5B give examples of source and receiver pairs (S, R) that produce reflections from a given location inside the body with a given reflection angle. Figure 5A shows a schematic cut-away perspective view 500 of a cylindrical body part, with an image point 502 inside the body. A source 504 on the surface of the body optionally transmits ultrasound waves into the body, optionally in all directions, and ray optics path 506 reaches point 502. Although the rays will be straight lines if the velocity model has a sound speed that is uniform everywhere inside the body, in general the velocity model will not have a uniform sound speed, and the rays will be curved, as shown in FIGS. 5A and 5B. A receiver transducer 510 receives ultrasound, optionally from all directions inside the body, and ray optics path 508 reaches transducer 510 from point 502. Path 506 and path 508 intersect at point 502 at an angle 2α_i, corresponding to a reflection angle α_i, so source 504 and receiver 510 constitute one of the source-receiver pairs that transmit ultrasound to point 502, and receive reflected ultrasound from point 502, with a reflection angle α_i . In general there are many source-receiver pairs that transmit ultrasound to point 502, and receive the reflected ultrasound with the same reflection angle α_i . Some of these pairs come from different imaging directions, where the imaging direction is defined at the direction of a line starting at reflection point 502 and oriented halfway between the direction of path 506 from the source and the direction of path 508 to the receiver. For example, source 512 transmits ultrasound along a ray path 514, that reflects from point 502 with a reflection angle α_i , and follows ray path 516 to receiver 518. The imaging direction for source-receiver pair 504 and 510 is almost opposite to the imaging direction for source-receiver pair 512 and 518. The two imaging directions are indicated by the solid arrows in view 500. There are also other source-receiver pairs that have the same imaging direction and the same reflection angle as source-receiver pair 504 and 510, for the same reflection point 502, but with the plane of the reflection oriented at a different angle.

There are also many source-receiver pairs that transmit ultrasound to point 502 and receive reflected ultrasound from point 502 with different reflection angles. For example, FIG. 5B shows a cutaway perspective view of the same cylindrical body part with the same point 502 inside the body. Source 522 transmits ultrasound along a ray path 524 to point 502, which reflects with a reflection angle 0C3, smaller than α_i, onto ray path 526 that is received by receiver 528. Source 530 transmits ultrasound along a ray path 532 to point 502, which reflects with a reflection angle o , larger than α_i , and travels along ray path 534 to receiver transducer 536. It should be noted that a given source-receiver pair (S, R) transmits and receives ultrasound from many different points inside the body at a given reflection angle. For example, if the velocity model has a uniform sound speed, so that the rays are straight lines, then the points for which a given source-receiver pair has a reflection angle of 45° will be the points, inside the body, on the surface of a sphere whose center is midway between the source and the receiver, and whose diameter is the distance between the source and the receiver.

If a plurality of closely spaced source transducers transmit ultrasound simultaneously with a well-defined phase relation between them, an option for some embodiments of the invention described above for 308 in FIG. 3, so that the transmitted beam is limited in angle, or if a single source transducer produces a transmitted beam of ultrasound that is limited in angle, then the only source-receiver pairs (S, R) that can transmit and receive a signal that reflects from a given point (r, θ, z) in the field of view, are the source-receiver pairs for which that point falls within the angular range of the transmitted beam of the source. In some embodiments of the invention with transmitted beams of limited angular width, if a given point (r, θ, z) falls within the angular range of the transmitted beam but near its edge, then the received signal for that source-receiver pair is included when finding M(r, θ, z, a) for that (r, θ, z), but with lower weight, for example in proportion to the angular profile of the transmitted beam. Alternatively, the received signals for all the sourcereceiver pairs are given the same weight when finding M(r, θ, z, a) for that (r, θ, z).

To find the migrated reflectivity M(r, θ, z, a) for a given point (r, θ, z) and reflection angle a, the received signals are summed for all the source-receiver pairs (S, R) for which a ray path from source S reaches point (r, θ, z) and reflects into another ray path at a reflection angle a, that reaches receiver R, each signal evaluated at a time TS,R(T, θ, z). In some embodiments of the invention, instead lumping together all of the signals for which the reflection angle is a, the signals are further sorted into groups according to one or both components of the angle of orientation of the plane defined locally by the ray paths of the transmitted wave and the reflected wave near point (r, θ, z), and the migrated reflectivity M is defined as a function both of a and of these one or two components of the orientation angle. That procedure may be useful especially if the body tissue at that (r, θ, z) has an anisotropic sound speed, which can be measured by the dependence of M on these components of the orientation angle.

Each signal in the sum is weighted by a factor fi equal to cosa, and by a factor f2 proportional to the volume of the grid element or voxel associated with the point (r, θ, z). Each signal optionally has the same envelope, which includes at least a few wave periods of the central ultrasound frequency, so that different signals in the sum can interfere with each other constructively or destructively, depending on their relative phases. The time TS,R(T, θ, z) at which each signal is evaluated is, for example, the time relative to a transmission time defined for that signal, which is optionally defined in the same way relative to the envelope for each signal, for example the transmission time is defined as a specified time in the middle of the envelope, where the amplitude of the transmitted signal is fairly high, for example about as high as it gets, or the transmission time is defined as a time about halfway between the starting time and end time of the signal, or the transmission time is defined as the starting time of the signal.

The reflectivity image at each point (r, θ, z), which is a sum of the migrated reflectivities M(r, θ, z, a) for the different values of the reflection angle a, is found by summing the recorded signal data from many different source-receiver pairs (S, R), each signal evaluated at a time TS,R(T, θ, z). This sum includes wave components that were not reflected from this point. But in the sum of signals contributing to the reflectivity at a given point (r, θ, z), the wave components that reflect from that point interfere constructively, adding up, while the other wave components interfere destructively, largely cancelling out.

At 406, the data in the angle gathers, the migrated reflectivity data M(r, θ, z, a), is optionally used to find a correction to the velocity model. An exemplary method for doing this will be described in FIG. 7A. Briefly, the method relies on the expectation that, if the velocity model is correct, then the reflectivity M(r, a) for a given (θ, z) should not have a systematic dependence on reflection angle a in a way that would be expected if there were an error in the velocity model. If the velocity model is correct, we expect that any local peak in M(r, a) as a function of r, caused by a feature at a distance beneath the surface, will occur at the same r for different values of a, or at least the value of r at which the peak in M(r, a) occurs should not have a steady drift with a. The main effect of an error in the velocity model is that the travel time along a path going to a point (r, θ, z) at an oblique angle a, relative to the travel time along a path going to the same location at a different angle, for example 0, will be calculated incorrectly, which will make M(r, a) systematically vary with a in a certain way. This can be used to correct the velocity model. Details of how this is done are given below in the descriptions of FIGS. 7 A, 11A, 1 IB, 12A and 12B.

At 408, an evaluation is made of whether the corrections to the velocity model are converged, so that the iterative process of correcting the velocity model, and using the new velocity model to find corrected angle gathers M(r, θ, z, a), can be stopped. The convergence criteria can include, for example, whether the latest corrections are sufficiently small, and whether the number of iterations so far is larger than a maximum number of iterations that will be allowed. If the velocity model is not converged, then the method returns to 404, and new angle gathers are generated, using the corrected velocity model to calculate the ray paths. If the velocity model is judged to be converged at 408, then pre-stack attributes are optionally found at 410. As used herein, “attributes” refer to characteristics of body tissue, including characteristics other than reflectivity, in the field of view as a function of location, that are potentially of medical interest, and that can be generated as images. Examples include the sound speed, the shear wave speed, the magnitude and orientation of sound speed anisotropy, Poisson’s ratio, and density. Here “characteristics of body tissue” includes quantities that are mathematically derived from quantities that represent direct physical characteristics of body tissue, even if the derived quantities do not have a simple physical meaning. Examples of mathematically derived quantities include time derivatives, gradients, spatial and temporal filtering, and complex traces of the quantities they are derived from. Pre-stack attributes are attributes that can be found directly from the angle gather data M(r, θ, z, a), using methods that will be described below. It should be noted that angle gathers can be useful in generating images of certain attributes of medical interest, such as Poisson’s ratio and S-wave velocity, even if they are not used to correct the velocity model. See FIG. 13 and the section “Pre-stack and Post-stack Attributes” below, for details on how pre-stack and post-stack attributes are found.

At 412, the main reflectivity image is found from the corrected angle gather data M(r, θ, z, a). With the velocity model corrected, M(r, θ, z, a) should not have a systematic travel time dependence on a of the type that would be produced by a error in the velocity model. By integrating the migrated reflectivity M(r, θ, z , a) over a, which is referred to as “stacking” the partial reflectivity images M(r, θ, z, a) over different values of a, a main reflectivity image M(r, θ, z) is found. This main reflectivity image M(r, θ, z) is similar to, but distinct from, the P-wave impedance reflectivity image, which is M(r, θ, z, a=0). The main reflectivity image may look similar to M(r, θ, z, a) for typical values of a, but with reduced noise. This image shows the main structure and texture of body tissues in the field of view, such as organs, bone, muscle, etc., which may be immediately recognizable by medical personnel viewing the image.

At 414, post-stack attribute images are found from pre-stack attributes, for example from the corresponding pre-stack attributes, for example using the inversion procedure described in detail at 1312 of FIG. 13. Post-stack attributes are also optionally derived by mathematical manipulation of other post-stack attributes, for example as described below in 1316 of FIG. 13.

At 416, the main reflectivity image M(r, θ, z) is outputted. At 418, additional images, of pre-stack and/or post-stack attributes, are optionally outputted. At 420, two or more images, for example the main reflectivity image and one or more attribute images, are optionally combined for multi-image visualization. For example, an attribute image, which might provide information about whether a tissue is cancerous, is superimposed on the main reflectivity image, which shows the tissues in a way that they can be immediately identified, so that someone viewing the image can immediately see where the cancer is located.

Although FIG. 4 describes a migration method in which angle gathers are used, there are other similar types of gathers that can be used instead, to provide similar information for correcting the velocity model. For example, instead of sorting the signal data by reflection angle a, as described in FIG. 4, the signal data can be sorted by the offset X, the distance between the source transducer and the receiver transducer. So the migrated reflectivity M(r, θ, z, X) can be used instead of M(r, θ, z, a) to correct the velocity model, and to find the main reflectivity image M(r, θ, z) by integrating over X. Using offset instead of reflection angle has the potential advantage that the offset between a given source transducer and receiver transducer is known independently of the velocity model and independent of the image location (r, θ, z), while the reflection angle for a given source and receiver and a given location (r, θ, z) has to be calculated using the velocity model and ray-tracing software. Using reflection angle a has the potential advantage that the reflectivity as a function of a can provide information on physical quantities such as Poisson’s ratio, S-wave velocity, and density. Other similar types of gathers, for example gathers that provide equivalent information for correcting the velocity model, can also be used.

Figure 6 schematically shows a perspective view 600 of an exemplary coordinate system (r, θ, z) used for generating an image of the field of view from ultrasound signal data, using the system of FIG. 2 and the methods of FIGS. 3 and 4. The coordinate system shown in FIG. 6 is an elliptic cylindrical coordinate system, with a straight z-axis, and with the outer surface of constant r shaped like an ellipse of the same shape and size and orientation at each value of z. Optionally, this coordinate system is chosen so that the z-axis corresponds approximately to a cylindrical axis of a roughly cylindrical body part that is being imaged, for example the trunk, and the elliptic crosssections correspond approximately to axial cross-sections of the body part, and the outer surface corresponds approximately to the surface of the body part.

It should be noted that in the coordinate system of FIG. 6, the lines of constant θ and constant z, for example lines 602, are straight lines, or nearly straight lines, that are normal or nearly normal to the outer surface of constant r, as well as to the other surfaces of constant r. The angle gathers, or whatever type of gathers are being used, follow paths of constant θ and constant z, so with this coordinate system the paths of constant θ and z will be straight lines, approximately normal to the surfaces of constant r, and approximately normal to the outer surface of the body, to the extent that the outer surface of the body roughly conforms to the outermost surface of constant r. Although the coordinate system of FIG. 6 is only one of many cylindrical coordinate systems that can be used for pre-stack depth Kirchhoff migration, there may be advantages to using a coordinate system for which the surface of the body is approximately a surface of constant r, and for which gathers extend into the body perpendicular to the surface. Another potential advantage of this coordinate system, which it shares with other cylindrical coordinate systems, is that the gathers get closer together in the 0 direction as they go deeper into the body, so smoothing over 0 can provide information about nearby gathers that potentially may compensate for the loss of accuracy that may occur further away from the surface where the transducers are located.

It should be noted that with this coordinate system, if the elliptical cross-sections are not circles, then the lines of constant θ and z do not all end at the z-axis, as they would for a circular cylindrical coordinate system, but they end at different locations along the major axis of the elliptical cross-section. Optionally, in this coordinate system, each of the surfaces of constant r is spaced at a constant distance inward from the outermost surface of constant r, with r increasing linearly with distance in from the outermost surface of constant r. Then r will be a measure of distance inward from the outermost surface of constant r, and will be approximately a measure of distance inward from the surface of the body part, to the extent that the outermost surface of constant r approximately matches the surface of the body part. However, since the body part will in general not be shaped exactly like an elliptical cylinder, in general the value of r at the surface of the body will vary with θ and z. In this coordinate system, r increases going inward from the outermost surface of constant r, and from the surface of the body. The value of r at the beginning of an angle gather at a given θ and z, which may be at the surface of the body, will be called r_min(θ, z), since it is the minimum value of r for that angle gather. The surfaces of constant r will, in this case, generally not be the same shape as the outermost surface of constant r, and may not be ellipses at all. It should be noted that an angle gather need not extend all the way to the outer surface of the body, for example if the field of view, or the region of interest to be imaged, does not extend all the way to the outer surface of the body.

Additionally or alternatively, even in a coordinate system where the gathers, the lines of constant θ and z, are not approximately normal to the surfaces of constant r and to the outer surface, optionally as you go along each gather in one direction you get further away from the nearest point on the outer surface, and r increases in this direction. In this sense, r is a measure of distance into the field of view, away from the outer surface. Optionally, on nearby gathers r increases at similar rates as a function of distance along the gather, and surfaces of constant r are fairly smooth and well-behaved. Points on neighboring gathers, with nearby values of θ and z and the same value of r, will be spatially close to each other, and smoothing of the velocity model over a surface of constant r makes sense physically and may make the velocity model more accurate, as will described below in the description of FIG. 7A. This may be true even if the rate of increase of r as a function of distance along a gather is very different for gathers that are not near each other.

Alternatively, the coordinate system is an elliptical cylindrical coordinate system, but all the surfaces of constant r are ellipses with the same eccentricity centered on the z-axis, and the lines of constant θ are straight lines going from the outermost surface of constant r to the z-axis. In such a coordinate system, the surfaces of constant r have a simpler analytic form, which is a potential advantage, but the lines of constant θ and z will generally not be normal to the surfaces of constant r, and r will not generally be a measure of distance from the outermost surface of constant r.

Alternatively, the (r, θ, z ) coordinate system need not have this form, but any three- dimensional coordinate system can be used, even including a Cartesian coordinate system, which might be particularly appropriate if the transducers are located only on one side of the field of view. The surfaces of constant θ and z will define the locations of the gathers used by the migration algorithm, and in principle, especially with Kirchhoff migration, the gathers can be located along any array of paths going into the field of view from the surface of the body part, and the migration algorithm will still work, as long as signal data is available from a complete enough array of transducers. But some coordinate systems may make the migration more accurate and/or more efficient, as discussed above.

In some embodiments of the invention, the coordinate system conforms to the shape of the body part surface, for example the body part surface is a surface of constant r. Or the outer surface of constant r at each value of z is an ellipse, but the size and/or eccentricity of the ellipse varies with z, to better match the surface of the body part. However, using a coordinate system with a simple analytic form, such as a circular or elliptical cylindrical coordinate system, has the potential advantages that calculations may be simpler, and that the coordinate system does not have to be adapted for each person being imaged, or to the posture of the subject during imaging.

For the case of an approximately cylindrical body part, especially with transducers located all the way around it in an azimuthal direction, for example, it may be convenient to use θ as a periodic coordinate in the azimuthal direction, to use r as a depth coordinate approximately measuring distance into the field of view from the surface of the body part, and to use z as an axial coordinate measuring distance along the axis of the body part.

Using angle gather data to correct the velocity model

Figure 7 A is a flowchart 700 showing a method of finding corrections to the velocity model, from the data in angle gathers based on ultrasound signals transmitted through a field of view of a body part, between different pairs of source transducers and receiver transducers mounted on the surface of the body part, according to an exemplary embodiment of the invention. Flowchart 700 provides more details of the loop that comprises 404, 406 and 408 of flowchart 400 in FIG. 4. In describing flowchart 700, the sound speed V of the velocity model will be considered a scalar, but the method can be generalized for the case where the sound speed may be an anisotropic tensor, for example in bone and in some soft tissue, and for the case where shear waves are considered. Optionally this is done, for example if the field of view includes tissue for which the sound speed is expected to be anisotropic, and/or tissue in which shear waves are expected to occur. Canning and Malkin, cited above, describe the case of anisotropic sound speed.

It should be understood that the velocity model V(r, θ, z), as well as other quantities that are calculated as functions of position (r, θ, z) in the field of view, such as angle gathers, are optionally defined only at discrete grid points (r, θ, z) in the (r, θ, z) coordinate system, and the values of these quantities at other locations (r, θ, z) may be determined if needed by some kind of interpolation method, such as linear, quadratic, or cubic spline interpolation. Equivalently, these quantities may be defined as analytic functions for all real values of r, θ, and z, for example using cubic spline coefficients.

At 702, the precise positions of the transducers are found. Knowing the positions of the transducers is generally needed for image reconstruction. For example, the procedures in the rest of flowchart 700, which include finding ray optic paths between each transducer and different points in the field of view, and finding travel times of ultrasound signals along these paths, assume precise knowledge of the positions of the transducers, in order to get accurate results. In particular, the positions of all the transducers, in r, θ, and z, are assumed to be known to within a fraction of a wavelength of the ultrasound being used. For example, if the ultrasound has a frequency of 1 MHz, then the wavelength is about 1.5 mm, and the positions of the transducers are assumed to be known to within about 0.75 mm, or even much better. This is difficult to accomplish, if the transducers are mounted on the surface of a soft body part, such as the trunk, whose shape may be constantly changing, for example due to breathing. The inventor has found that it is possible to determine the positions of the transducers to such a high precision, by a machine learning algorithm, using the recorded ultrasound signals as input. Optionally, sensors, such as strain sensors and/or angle sensors, in the flexible stretchable frame of the array of the transducers, provide an initial estimate of the positions of the transducers, and this estimate is used as supplementary input to the machine learning algorithm.

Figure 9A describes an exemplary convolutional neural network (CNN) that can be used as the machine learning algorithm, and FIGS. 10A and 10B describe how the CNN is trained. The CNN, already trained, or another kind of trained machine learning algorithm, is run at 702, using as input the same recorded signal data, or at least a portion of it, that will be used later in flowchart 700 to find ray optics paths and travel times in order to find the angle gathers. The CNN will then find precise positions of all the transducers, which may be used, among other purposes, in order to find the ray optics paths and travel times to adequate precision. If the signal data is recorded over an extended period of time, during which the positions of the transducers may change significantly, then the CNN optionally finds the position of each transducer at each of the times the different signals that involve that transducer are recorded. Using the CNN to find the positions of the transducers may involve making multiple applications of the CNN, if each application of the CNN only finds the positions of a subset of the transducers, as described below, and then using software to integrate and reconcile the outputs of the different applications of the CNN, including taking into account possible changes in positions of the transducers over time. The CNN only has to find precise positions of the transducers once, at the beginning of flowchart 700, because for the rest of flowchart 700 the same recorded signal data is used. Only if new signal data is recorded, is there any reason to use the CNN again to find the positions of the transducers.

Optionally, for example if there is reason to be concerned that the positions of the transducers may have changed during the time when the signals were transmitted and recorded, the CNN is used to determine whether and how the positions of the transducers have changed over this time. For example, the CNN can be used several times, each time using only signal data recorded during a relatively small fraction of the total time that signal data is recorded, during which the transducer positions are not expected to change significantly. Precise positions of the transducers can then be determined, for example by interpolation, at the time when each signal was transmitted, and these changing positions can be taken into account when using the recorded signal data to reconstruct an image. This might be done, for example, when imaging the trunk, if the time needed to record all of the signals is too long for the subject to hold their breath, or if it is difficult for the subject to hold their breath or to remain still enough so that the transducer positions do not change significantly, for example by more than a fraction of 1 mm.

Box 704 lists the elements of a procedure for performing migration of the recorded signals us,R(t), in order to obtain information about each point (r, θ, z) in the field of view. Each us,R(t) is the recorded signal, for example the perturbed pressure, as a function of time t, transmitted by a source transducer S and received by a receiver transducer R. Here, the time t is measured relative to a transmission time of the signal by the source transducer, defined in the same way for all signals. For example, the transmission time is defined as the start time of the signal, or as a time halfway between the start and end time of the signal, or as a specific time in the envelope of the signal that has relatively high amplitude, for example about as high as the amplitude gets. Optionally, to ensure that the time t in the recorded signal data us,R(t) is defined relative to the transmission time of the signal by the source transducer, controller 212, in FIG. 2, synchronizes the transmission of the signal by the source transducer and the recording of the received signals by the receiver transducer, for example by sending a command to the source transducer S to start transmitting the signal, and at the same time sending a command to each receiver transducer to start recording the received signal, with t = θ at the start of recording the signal. In that case, the transmission time would be defined as the beginning of the transmission signal. In order to define the transmission time at a different time, for example at a time At after the start of the transmitted signal, for example in the middle of the envelope of the transmitted signal, the receiver transducer instead optionally sets t = -At at the time it begins recording the received signal, or the receiver transducer optionally only starts recording the received signal at a time At after the start of the transmission by the source transducer, and sets that time as t = 0, although in that case data for some of the received signal may be lost.

The absolute time, for example the clock time, of transmission of the signal of each source transducer, however it is defined, is optionally recorded and saved as a “time stamp” associated with the signal transmission by that source transducer. As described below, this time stamp may be useful for finding the positions of the transducers as a function of time by the CNN, or by software that integrates the output of different applications of the CNN. Since the time stamp is only going to be used for the purpose of finding the position of the transducers as a function of time, it does not need to have very high precision, but need only be precise to a time over which the positions of the transducers are expected to change significantly, for example by a small fraction of an ultrasound wavelength. This time might be on the order of a few milliseconds, if the surface of the body is moving a few centimeters per second, for example due to breathing. By contrast, the synchronization of the source transducer transmitting the signal, and the receiver transducer recording the signal, should be precise to within a small fraction of an ultrasound wave period, for example to within 0.1 microseconds.

More precisely, migration is used to obtain information about each of a set of discrete points in the field of view that is being used for the image, for example points located on a discrete grid, corresponding to voxels of the image. The type of migration described in flowchart 700 is pre-stack depth Kirchhoff migration, and it is particularly well suited for a field of view with an irregularly shaped outer surface on which transducers are located, such as the field of view of a body part that is being imaged. It is also possible to adapt Kirchhoff migration to different coordinate systems, such as a cylindrical coordinate system that might be convenient to use for a cylindrical body part. And Kirchhoff migration may be especially suitable for accurately taking into account the precise positions of the transducers found by the CNN, and their small changes over time. However, other types of migration can also be used.

The pre-stack depth Kirchhoff migration procedure described in box 704 generates an angle gather at each of a discrete set of pairs (θ, z). Each angle gather can be thought of as starting at a point with those values of θ and z at the surface of the body part, and extending along a path of constant θ and z, going into the field of view. The coordinate r increases going deeper into the body along a path of constant θ and z. For a given angle gather defined by a pair (θ, z), r ranges from its minimum value r_min where the angle gather begins, for example at the surface of the body part, to a maximum value r_max inside the body where the angle gather ends. For example, if the coordinate system is a cylindrical coordinate system, then the angle gathers may end at the z-axis. The r coordinate need not be constant on the surface of the body part, so r_min may be different for different angle gathers ( θ, z). Some of the angle gathers may extend deeper into the body than other angle gathers, and r_max may also be different for different angle gathers (θ, z). The data in each angle gather is a function of r and of a reflection angle a. The angle gather data is used, in the second part of flowchart 700, to obtain a correction to the velocity model.

It should be noted that other types of gathers, using a variable other than reflection angle, can be constructed with Kirchhoff migration instead of angle gathers to obtain the correction to the velocity model. For example offset gathers can be used, where the variable is the offset, or distance, between a source transducer and a receiver transducer, and the procedure of flowchart 700 can be modified, in a straightforward way, to use offset gathers or other similar types of gathers instead of angle gathers.

The procedure in box 704 is performed for each of the points (r, θ, z) in the field of view that are being used for the migration. For each angle gather ( θ, z), the procedure is performed for each r in the range of that angle gather, resulting in a migrated reflectivity M(r, θ, z, a). This migrated reflectivity is then used for the procedure in box 716, to find a correction AV to the average velocity Vims for each r in each angle gather. The corrected average velocities V_rms for all the angle gathers are then used in 736 to find a corrected interval velocity model V(r, θ, z), and the corrected velocity model V(r, θ, z) is then used, in an iteration loop, to perform all of the procedures in boxes 704 and 716 over again. Within each iteration of this iteration loop, there is a potential advantage in performing all of the procedure of box 704 for a given angle gather (θ, z), and then to perform all of the procedure of the following box 716, for the same angle gather (θ, z), and then going back to the beginning of box 704 for the next angle gather (θ, z), because this greatly reduces the storage requirements of the computation. Alternatively, within each iteration, each element of the procedure in boxes 704 and 716, or some part of the procedure, is performed for all the angle gathers, storing the results, before going on to the next element of the procedure.

At 706, for each of the points (r, θ, z) in the field of view, paths of ultrasound waves are calculated, using ray optics, from point (r, θ, z) to each transducer on the surface of the body part, using the current velocity model V(r', θ', z'). The sound speed V has been written here as a function or r', θ', and z' , rather than as a function of r, θ, and z, to show that the ray optics path from a given point (r, θ, z) to a given transducer depends on the sound speed V at other points in the field of view, not just the sound speed at that point. The velocity model affects the ray optics paths of the ultrasound, due to refraction. For example, if the velocity model has a uniform sound speed V, independent of r, θ, and z, then the ray optics paths will all be straight lines. Initially, before the velocity model has been corrected, an initial velocity model is used. Optionally, the initial velocity model has a uniform sound speed V, for example equal to 1540 m/s, the sound speed in water at body temperature and at the salinity of blood, which is a fairly good approximation to the sound speed in soft body tissue. This choice may be especially useful if the field of view is mostly soft tissue, as it is in the trunk, and in the fleshy parts of the arms and legs. Particularly for a field of view that includes bone as a significant fraction of the volume, an alternative initial velocity model is optionally used, with a higher sound speed, appropriate for bone, in parts of the field of view that are expected to be bone. Such an initial velocity model is also optionally used even in the trunk, or in the arms and legs, modeling which parts of the field of view are thought in advance to be bone.

Also at 706, the travel time for an ultrasound wave along each path, from the point (r, θ, z) to each of the transducers, is calculated, optionally in the ray optics limit, for example by integrating 1/V(r', θ', z') along the path. This calculation also depends on the velocity model V(r', θ', z'), and the initial velocity model is used the first time, before it has been corrected.

At 708, for each source-receiver pair (S, R) for which signal data us,R(t) has been recorded, and for each point (r, θ, z) in the field of view, the reflection angle αS,R(r, θ, z) at the point is calculated. Here S and R are identifiers, for example numbers, labeling the source transducer and the receiver transducer. Each transducer in the array has a unique identifier. The reflection angle ocs,R(r, θ, z) is half of the angle between the ray optics path from point (r, θ, z) to the source transducer S, and the ray optics path from point (r, θ, z) to the receiver transducer R.

Also at 708, the total travel time TS,R(T, θ, Z), from the source transducer S, to the point (r, θ, z), and from the point to the receiver transducer R, is found, by adding the travel time from the point to the source transducer, and the travel time from the source to the receiver transducer, that were found at 706. At 710, a discrete set of reflection angles { α_i } is specified, for example between θ and a maximum reflection angle α_max that will be used for the angle gathers. Optionally this is done only once, for example when the software is written, and the same set { α_i } is always used, or one or more input parameters specify the set { a/} . In seismic acoustic imaging using Kirchhoff migration with angle gathers, the maximum reflection angle α_max used is typically between 30° and 50°. It is difficult to use reflection angles much greater than that for seismic imaging, for most underground depths of interest, because of the limited lateral extent of the array of acoustic sources and geophones on the surface of the Earth, and because typically, in seismic imaging, at larger reflection angles there is a large amount of interference from “Head Waves.”

Figure 7B shows a schematic cross-sectional view 750 of an image point 752 in a field of view 754 used for seismic imaging. In order to obtain a reflection angle of 80°, an acoustic source 756 and a receiver 758, on a flat surface 760 of the Earth, would have to be located far from the field of view, compared to the dimensions of the field of view. Ray path 762 from source 756 to image point 752 makes an angle of 160° with a ray path 764 from image point 752 to receiver 758. In a schematic cross-sectional view 766 of an image point 768 in a round body part 770, with a field of view for ultrasound medical imaging that corresponds to the entire cross-section of the body part, an ultrasound source 772 and an ultrasound receiver 774 are located on surface 776 of the body part, at locations for which the reflection angle is 80°. A ray path 778 from source 772 to image point 768 makes an angle of 160° with a ray path 780 from image point 768 to receiver 774. Because of the curvature of body surface 776, which completely surrounds field of view 770, the source and receiver can be located relatively close together, compared to the dimensions of the field of view.

Potentially, for medical ultrasound imaging with transducers located all the way around the body part being imaged, α_max could be as high as 90°. A reflection angle A that is greater than 90° would be equivalent to a reflection angle of 180° - A which would be less than 90°, for an angle gather on the opposite side of the field of view. However, for reflection angles close to 90°, it may be difficult to get good data, because the path from the source to the point to the receiver will be almost the same as a direct path from the source to the receiver, and the signal may be dominated by the direct transmission from source to receiver with no reflection, which may be much stronger than the part of the signal due to reflection from points in the field of view. For this reason, α_max is optionally chosen to be well below 90°, for example between 50° and 70°. Alternatively, α_max is chosen to be closer to 90°, for example between 70° and 80°, or between 80° and 85°, or between 85° and 90°. Optionally, α_max is greater than 90°, even up to 180°, so reflection angles a greater than 90° are included in the angle gather, instead of being assigned to reflection angle 180° - a for another angle gather on the opposite side of the field of view, or in addition to being assigned to the other angle gather. The discrete angles α_i are the reflection angles a at which the angle gather data M(r, θ, z, a) will be calculated, and are optionally a few degrees apart, for example 0.5° apart, or 1° apart, or 2° apart, or 5° apart, and the set has, for example, 10, 20, 50 or 100 discrete angles α_i. If the discrete angles α_i are too close together, so there are more of them for a given α_max, then each M(r, θ, z, a) will have contributions from fewer recorded signals, and M(r, θ, z, a) may have a high noise level. If the discrete angles α_i are too far apart, then it may be difficult to determine how M(r, θ, z, a) varies with a, because M(r, θ, z, a) may be evaluated at too few values of a for a given point (r, θ, z).

At 712, for each point (r, θ, z), for each source-receiver pair (S, R) for which signal data was recorded, it is determined which of the discrete angles α_i is closest to the reflection angle ocs,R(r, θ, z) found at 708.

At 714, for each discrete reflection angle α_i in the set, for each point (r, θ, z), the time derivatives of the recorded signals us,R(t), evaluated at t = TS,R(r, θ, z), are added together, for those pairs (S, R) for which it was determined at 712 that αS,R(r, θ, z) is closest to α_i. Each term in the sum is weighted by the factors fi and f2, where fi = cos( αS,R(r, θ, z)), which is optionally approximated by cos(α_i), and f2 is proportional to the volume of the grid element or voxel at (r, θ, z). The sum is the migrated reflectivity M(r, θ, z, αi).

It should be understood that the procedures in box 704, which are done for each point (r, θ, z) in the field of view, can be done in any order. For example, all of the procedures in box 704 can be done first for one point (r, θ, z), then for the next point, and so on, because what is done at each procedure of box 704 depends only on the results of calculations made previously for that point (r, θ, z), not on the results of calculations for other points. Doing all these calculations for one point (r, θ, z) at a time has the potential advantage of reducing memory requirements for the computation, because, for example, the data on the ray optics paths from a given point (r, θ, z) are optionally discarded once M(r, θ, z , α_i) has been calculated for that (r, θ, z) for all the α_i’s, Similarly, procedures 708 and 712, and optionally 714, may be done for only one source-receiver pair (S, R) at a time, and optionally some or all of these procedures are done for each pair (S, R), before doing them for the next (S, R). Alternatively, some or all of the procedures in box 704 are done for all points (r, θ, z), and/or for all pairs (S, R), before going on to the next procedure.

Box 716 lists the procedure for using M(r, θ, z, α_i , at each point (r, θ, z) in the field of view and each reflection angle α_i, to find a correction to the local average velocity V_rms(r, θ, z), which will be defined. This correction is found independently for each angle gather (θ, z). Once the corrected V_rms(r, θ, z) is found for all the points in the field of view, it is used to find a correction to the velocity model V(r, θ, z), as will be described below.

At 718, the vertical travel time t(r, θ, z) is calculated, for a sound wave to go from the surface of the body part at a given θ and z, down to the point (r, θ, z) at the same θ and z inside the body. It should be understood that use herein of the term “vertical travel time,” which is borrowed from the field of seismic imaging, and expressions such as “down into the body,” “beneath the surface of the body,” “top of the angle gather” and “bottom of the angle gather,” do not imply that the path from the surface of the body to a point inside the body, keeping the same θ and z, is literally vertical. If, as in the coordinate system shown in FIG. 6, the lines of constant θ and z are normal to the surfaces of constant r, and the outermost surface of constant r is approximately the surface of the body, then point (r, θ, z) will be approximately straight down into the body along the angle gather at that (θ, z) from the surface of the body. The vertical travel time is given by t(r, θ, z) = Sometimes this will be written herein as t(r), with the dependence on 0

and z understood.

At 720, the average velocity V_rms(r, θ, z) is found by taking the rms average of V(r, θ, z) along the path of constant (θ, z), from r_min, where the angle gather begins at the surface of the body part, down to r, but averaging over time t(r) instead of over r. So V_rms(r, θ, z) = . This may also be expressed as an average over r, V_rms(r, θ, z) =

At 722, the angle gather data Me,_z(t(r), a), for the angle gather at a given (θ, z), is defined as M(r, θ, z, a).

The rest of the procedures in box 716 constitute an iterative loop for finding a correction to the average velocity V_rms(r, θ, z) from the angle gather data Me,_z(t(r), a). This calculation is local to each angle gather (θ, z), and is done for each of the angle gathers (θ, z). Optionally, the loop is performed for one angle gather (θ, z) until the correction converges, and then the loop is performed for the next angle gather (θ, z), and so on, until the converged correction has been found for all the angle gathers. Alternatively, each element of the loop is performed for some or all the angle gathers, before going on the next element, or a group of the elements of the loop is performed for some or all the angle gathers, before going on the next element, or before going on to the next iteration. A potential advantage of iterating the whole loop for each angle gather until the velocity correction converges, before going on to the next angle gather, is that it may reduce memory requirements. Optionally, elements 718, 720, and/or 722, are also performed for each angle gather (θ, z), and the loop is iterated until it converges for that (θ, z), before elements 718, 720, and 722 are performed and the loop is iterated for the next angle gather (θ, z), which also has the potential advantage that it may reduce memory requirements.

At 724, 1 + Gsin²a is fitted to Me,_z(t, a) for each t in the gather. It should be understood that the best fits I and G are functions of θ, z, and t, so they will be written as Ie,_z(t) and Ge,_z(t). Physically, Me,_z(t, a) represents the reflectivity of a point in the field of view, at this (θ, z), that is located far enough below the surface of the body so that a sound wave would take a time t to reach it going straight down along the path of constant (θ, z). The reflectivity can be measured by a source and receiver centered at equal lateral distances in opposite directions from the top of the angle gather, the point at that (θ, z) on the surface of the body. The different lateral distances correspond to different reflection angles a, with the lateral distance from the top of the angle gather to the source or the receiver being roughly tana times the vertical distance, which is approximately Vrmst, if the surface of the body is treated as flat and horizontal, if the oblique ray paths of the sound waves are treated as straight lines, and if the angle gather, the line of constant (θ, z), is treated as vertical. If the velocity model is correct, and in particular if Vims is correct, then we expect Me,_z(t, a) not to have a large systematic drift with a, and G would then be small.

G is a measure of how much Me,_z(t, a) systematically varies with a, in the way that would be expected if the velocity model is wrong on average, in the vicinity of that angle gather. G can be used to find a local error in the velocity model, specifically an error in the average velocity V_rms. Although the local error in Vrms, calculated in this way, may not be very accurate, because the body surface is not flat and horizontal, the sound waves may not follow straight ray paths, and the angle gather may not be vertical, the sign of the error will usually still be correct, especially if the error is large, and an iterative procedure can usually be used to zero in on a correct velocity model, or to get very close to a correct velocity model.

At 726, the complex trace I of Ie,_z(t) is found, and the complex trace G of Ge,_z(t) is found. The complex trace of a function is defined as the function, plus i times its Hilbert transform. For example, the complex trace of cos(ot) is exp(iot). It should be understood that I and G are also functions of θ, z, and t.

At 728, a estimated local correction AV(r, θ, z) to the average velocity Vrms(r, θ, z) is found, from Ie,_z(t) and Ge,_z(t) and their complex traces:

where fo is the central frequency of the ultrasound. This equation is derived, for example, by Swan, cited above, and by Canning and Malkin, cited above, for the case of seismic acoustic imaging. An explanation for this result will be given below in FIGS. 11A, 1 IB, 12 A, and 12B.

At 730, Vnns(r, θ, z) is updated by adding AV(r, θ, z) to it. Alternatively, instead of updating Vrms at every cycle of the iteration, Vims is updated only at the end of all the cycles, when the correction is converged, but adding to Vrms all the corrections AV generated by all the cycles of the iteration. Although the result will not be exactly the same, because AV depends on Vrms, the result will be very close, because the total correction to Vrms will generally be only a small fraction of Vrms.

At 732, Me,_z(t, a) is corrected, replacing it by Me,_z(t + At, a), where At is the error in the implied vertical travel time t that would result from an error AV in Vrms, for a source and receiver that have an offset from the top of the angle gather that would result in a reflection angle a. The error in vertical time is

Δt = 2r sin²a AV/Vrms²

This expression for At is derived in the description of FIG. 11 A, below.

At 734, the loop begins at 724, using the corrected Me,_z(t + At, a). The loop iterations continue until a convergence criterion is satisfied, calculating AV and At and correcting Me,_z(t, a) each cycle, for all the values of a in the set {a/} and for all the values of r in that angle gather. The convergence criterion is, for example, that the correction AV falls below a minimum value for all values of r in that angle gather, for example that AV is less than a specified fraction of Vrms for all values of r, for example less than 0.5% or less than 1% or less than 2% or less than 5% or less than 10% of Vrms, or less than a smaller, a larger, or an intermediate fraction of Vrms. Optionally, the convergence criterion is 10% of the expected range of sound speeds in soft tissue, in the body part being imaged. For example, the sound speed in soft tissue of the trunk varies from 1440 m/s in subcutaneous adipose tissue to 1640 m/s for intervertebral disks, a range of 200 m/s. Sound speed in tendons is about 1750 m/s. Optionally, there is a maximum number of iterations allowed, for example 10 iterations, 20 iterations, 50 iterations, or 100 iterations, or a smaller, larger or intermediate number of iterations, and the iterations end at that number, even if AV is not small enough to meet the convergence criterion for all r in the angle gather.

After the loop between 724 and 734 has converged for all r, for all angle gathers (θ, z) in the field of view, the result will be a corrected average velocity V_rms(r, θ, z) for all points within the field of view. This corrected Vrms has been found based on the assumption that the apparent a- dependence of M(r, θ, z, a) for each angle gather (θ, z) is due entirely to errors in the average velocity model locally around that angle gather. This is indeed likely to be the dominant cause of the apparent a-dependence of M(r, θ, z, a), but errors in M(r, θ, z, a) can also be caused by neighboring angle gathers having different average velocity, since the ray paths for obliquely travelling sound waves pass through neighboring values of (θ, z), so the local corrected error Vims that results after iterating the loop between 724 and 734 may not be correct. However, when the corrected average velocity V_rms(r, θ, z) is used to find a corrected velocity model V(r, θ, z), with smoothing over neighboring values of θ and z, and when the procedure from 706 to 736 is iterated several times using the latest corrected velocity model V(r, θ, z), as will be described below, then the final corrected velocity model V(r, θ, z), which takes into account the θ and z dependence of the local Vims, is likely to be more accurate.

At 736, the average velocity V_rms(r, θ, z), which is an average over each angle gather from rmin, the top of the angle gather, down to r, is used to find the interval velocity V(r, θ, z), the velocity model at each point in the field of view. At 738, the new velocity model V(r, θ, z) is used to perform all the procedures from 706 to 736 over again, and this large iteration loop continues until a convergence criterion is satisfied for V(r, θ, z). For example, the largest change in V(r, θ, z) has to be smaller than a specified value, for example a specified fraction of that value of V(r, θ, z), or a specified fraction of an average of V(r, θ, z) over the whole field of view, or a specified fraction of the difference between the maximum and minimum V(r, θ, z) in the field of view. The specified fraction is, for example, 0.2%, or 0.5%, or 1%, or 2%, or 5%. Optionally, a maximum number of iterations is allowed in this large loop, regardless of how much V(r, θ, z) changes from one iteration to the next. In practice the inventor has found that V(r, θ, z) changes very little after between 2 and 5 iterations, and optionally a maximum of 3, 4 or 5 iterations is allowed. Alternatively, a greater maximum number of iterations is used, for example between 6 and 10, or greater than 10, or there is no maximum number of iterations, or a very large maximum number of iterations is used, designed only to indicate an error in the software if it seems to be stuck in an infinite loop.

In some embodiments of the invention, a procedure called “forced flattening” is used, after the iterations from 706 to 738 have been completed. The inventor has found that using this forced flattening procedure after the last iteration can result in an improved image quality. In forced flattening, Me,_z(t, a) for each angle gather (θ, z) is replaced by Me,_z(t + r_p(t,a), a), where r_p(t,a) is a time shift that is chosen to maximize the correlation between Me,_z(t + r_p(t,a), a) and me,_z(t), where me,_z(t) is a model trace equal to an average or typical Me,_z(t, a) for small reflection angles a, much less than 1 radian. To find r_p(t,a), the cross -correlation c(t, T) between me,_z(t) and Me,_z(t, a) is calculated for each reflection angle a, more precisely for each bin α_i of reflection angles that was used to create the angle gather data. If Me,_z(t, a), and me,_z(t) are treated as functions of t for all real values of t within the time limits of the signals, for example using cubic splines or some other basis functions with coefficients depending on the sampled values of the signals, then the crosscorrelation is given by

Then the time shift is the value of T for which is maximized. Here W(t') is a time

window, typically a time window for which for example a time window that is

symmetric around t'= 0, that is wide enough to sample at least a few ultrasound wave periods of Me,_z(t, a) and me,_z(t), or at least a few times the resolution time of the signals if that is greater than a wave period, for example at least 2 or at least 3 or at least 5 wave periods or resolution times, so that c(t, T) will be a meaningful measure of their local cross-correlation. But it is potentially advantageous not to make the time window too wide, for example more than 5 or more than 10 or more than 20 wave periods or resolution times wide, because that will smooth c(t, T) as a function of t, which will generally smooth r_p(t,a) as a function of t, which may make T_P(t,a) less accurate, and may make the resulting image less accurate. If Me,_z(t, a) and me,_z(t) are only defined on a set of discrete time values b, for example at a set of discrete times at which the signal data is sampled, then the cross-correlation c(t, T) may be defined as a sum over the discrete times, rather than as an integral over t', and c(t, T) will only be defined for those values of /, and for certain values of r for a given /. For example, if the discrete times b are equally spaced at intervals of a signal sampling time 6t, then b = to + i5t, and c(t, T) will be defined at T equal to integer multiples of 6t. Then c(t, T) may be given as

Here Wk is a time window, for example satisfying for example with Wk = W-k, and the

same considerations apply to the width of this time window as described above for the width of the time window W(l '). The sampling time 6t is, for example, at least a few times shorter than the wave period or the resolution time of the signal, for example at least 5 times shorter or at least 10 times shorter, to achieve a desired time resolution of the signal and a desired spatial resolution of the image. Then the time window will be much wider than 6t, and at least a few times wider than a wave period or a signal resolution time, and c(t, T) will be a meaningful measure of the local crosscorrelation between Me,_z(t, a) and me,_z(t).

When the ray optics paths are calculated with the revised velocity model V(r, θ, z), then the variation in V(r, θ, z) between different angle gathers (θ, z) is fully taken into account, in calculating Me,_z(t, a) for a given angle gather, even if the initial velocity model was uniform spatially. When V(r, θ, z) has converged, then the absence of a systematic a dependence in the Me,_z(t, a)’s, of the type expected from an error in the velocity model, may indicate that V(r, θ, z) is indeed accurate throughout the field of view, and that the absence of such a systematic a dependence is not just a result of the approximations made (for example straight line oblique ray paths, flat horizontal body surface, vertical angle gather) in the expressions for AV and At.

Optionally, when the interval velocity V(r, θ, z) is found from the average velocity V_rms(r, θ, z), V(r, θ, z) is smoothed over θ and/or z. This has the potential advantage of making V(r, θ, z) more accurate in the field of view as a whole, because the value of V(r, θ, z) calculated for a given point will depend not just on V_rms(r, θ, z) at nearby values of r for the same angle gather (θ, z), but also on Vmis(r, θ, z) for nearby angle gathers. This improves accuracy, because we expect V(r, θ, z) to be relatively uniform over small volumes that comprise the same type of tissue. As we go deeper into the body along the path of constant (θ, z) associated with an angle gather, the calculation of V(r, θ, z) from values of V_rms(r, θ, z) in the same angle gather will tend to have greater and greater errors, because V_rms(r, θ, z) may be less accurate at greater values of r, for example because it is based on data of ultrasound traversing longer paths through the body. The errors on neighboring angle gathers at the same r may tend to be uncorrelated with each other, and by averaging over nearby angle gathers at the same r, the errors may be reduced. This is especially true of smoothing over θ, because the paths of constant (θ, z) for angle gathers at neighboring values of 0 tend to get closer together as they go deeper in the body, as may be seen in FIG. 6. As we go deeper into the body, and the errors in V_rms(r, θ, z) on each angle gather may increase, the nearby angle gathers get closer together, and may provide more reliable information about the correct value of V_rms(r, θ, z) on that angle gather. By averaging over z and especially over θ at each value of r, the errors in V(r, θ, z) may be kept small, instead of increasing at greater r.

In an exemplary method of finding V(r, θ, z) from V_rms(r, θ, z), at the top of each angle gather, where r = r_min(θ, z), V(r, θ, z) is set equal to V_rms(r, θ, z). That may be done because V_rms(r, θ, z) is an average of V(r', θ, z) between r_min and r, along the path of constant (θ, z) for that angle gather. So at r = r_min, only the value of V at r_min contributes to the average, so the average is equal to V. The values of r at successive grid points along the path of constant (θ, z) for an angle gather, going into the body from the surface, will be designated n, r2, n, ..., up to r_max, the value of r at the far end of the angle gather inside the body, and r_min will be designated ro. It should be understood that n for a given value of the index i may refer to different values of r for different angle gathers, because r may have different values at the surface of the body for different angle gathers. In some embodiments of the invention, the top of an angle gather is always set at a grid point θ, z),

even if the outer surface of the body, at that (θ, z), does not correspond to a grid point. In that case the angle gather for that (θ, z) may begin slightly inside the body, at the first grid point for that (θ, z) that is located inside the body. This may mean that, for that (θ, z), the reconstructed image does not extend all the way to the surface of the body but only to the outermost grid point for that (θ, z) that is still inside the body. This may not matter very much if that grid point is near the surface of the body. In some cases, particularly if the surface of the body nearly corresponds to a surface of constant r, r_min may have the same value for all gathers (θ, z), and it may be convenient to define this value of r_min as r = 0, so that r is at least approximately a measure of distance into the body from the surface. In other cases, especially if the surface of the body is very irregular, so that it is not convenient to define a coordinate system where the surface of the body is at nearly the same value of r for all angle gathers (θ, z), r_min may have different values, possibly very different values, for different angle gathers (θ, z), and consequently, n for i = 1, 2,... may also have different values for different angle gathers (θ, z). For a given angle gather (θ, z), the interval velocity V(n) at each n is given by the following formula, given by C. H. Dix, C. H., “Seismic velocities from surface measurements,” Geophysics 20, 68-86 (1955).

It should be understood that in this formula, V, Vims and t are all functions of θ and z, as well as r, but since this formula is used only for a given angle gather, with a given value of θ and z, the 0’s and z’s have been suppressed in this formula.

Optionally, after the formula above has been used to calculate the V(r0’s for each of the angle gathers, the V’s are smoothed over θ and/or z, for at least some values of r in the field of view, for example for those values of r that are included in all or most of the angle gathers, or for all values of r in the field of view. Optionally, the values of V(r, θ, z) are calculated for a given value of r, for all the angle gathers, using the formula above, and the smoothing of V over θ and/or z is done for that value of r, before calculating V(r, θ, z) for the next value of r. Alternatively, the values of V(r, θ, z) are calculated using the formula above, for all the points in the field of view, and V(r, θ, z) for each value of r is then smoothed over θ and/or z. Optionally, V(r, θ, z) is also smoothed over r. The smoothing may be done by any known smoothing or filtering algorithm.

Convolutional neural network for finding positions of transducers The convolutional neural network (CNN) that finds precise positions of the transducers uses as input the recorded signals for many source-receiver pairs, using the timing of received signals transmitted directly from the source to the receiver, relative to the time of transmission, for example the time of transmission of the very beginning of the signal or the time of transmission of the exact middle of the signal, and/or the timing of received signals transmitted indirectly after reflecting from body tissue, relative to the time of transmission. However, for the large number of transducers, and recorded signals, needed to generate a high resolution image of a large field of view of a body part, for example an image with a resolution better than 2 mm, with a field of view covering several major organs in the trunk, it may not be feasible, with current computer technology, to use all of the recorded signals at once as input for the CNN. Furthermore, it may not be necessary to use all of the recorded signals at once, or to use all or even most of the recorded signals, as input for the CNN, in order to get accurate results for the positions of all the transducers. Optionally, a subset of the source transducers and receiver transducers are selected to provide input for the CNN, and signals from all the combinations of those source transducers and receiver transducers are used as input to the CNN.

Because it is generally not feasible to run the CNN with more than a small fraction of all the source transducers and receiver transducers in the array, the CNN may be run many times, each time using signal data from a different subset of the source transducers and receiver transducers in the array. For example, each subset consists of transducers located in a different patch of the array, small compared to the overall dimensions of the array, and with the patches overlapping, and with each subset optionally including all or most of the transducers in the area covered by that patch. The output of the CNN using data from each of the different overlapping patches, which gives the relative positions of the transducers in each patch, is combined to give the relative positions of all the transducers in the array. Alternatively, the combined output of the applications of the CNN gives the relative positions of a large number of widely distributed transducers in the array, not too distant from each other, consisting of only a small fraction of all the transducers in the array, and the positions of other transducers located between them can be estimated by interpolation. Optionally, instead of the different subsets each extending over a relatively small patch of the array, the different subsets each consist of a relatively broad area of the array, with transducers distributed sparsely over the area of the array that they cover, and with many different areas overlapping a given location in the array. That arrangement has the potential advantage that each subset includes transducers that are relatively far apart, which potentially provides good input to the CNN, but the number of transducers in each subset is small enough so that it is feasible to run the CNN. Optionally both types of subsets of transducers, and optionally also intermediate cases, are used to train the CNN, which may make the CNN work well in the application phase using either type of subset of transducers.

Depending on which signals are used as input to the CNN, the input signals may include a contribution from direct transmission of ultrasound from the source to the receiver, without reflection of ultrasound from body tissue, and may include a contribution from ultrasound that reflects from body tissue in going from the source to the receiver. The inventor believes that it is potentially useful to use input signals that have substantial contributions both from direct transmission of ultrasound, and from ultrasound that reflects from body tissue, for determining the relative positions of the transducers.

The CNN, being a program that can be trained by machine learning, is initially run by the programmers in a training phase, using a wide variety of training examples that resemble the input data that it will be exposed to when it is used by the end users. The training phase, which is described in detail in FIGS. 10A and 10B, is used to optimize a large number of free parameters, referred to as weights, so that it will learn to perform the task it will be used for, using recorded ultrasound signals from a subject to determine the precise positions of the transducers.

Once the CNN has been adequately trained, by one or more specified criteria as described in FIGS. 10A and 10B, the CNN can be used by end users in an application phase, using as input the recorded ultrasound signals of subjects who are undergoing ultrasound imaging, using the system and methods described above in FIGS. 2, 3, 4 and 7A. Typically, in the application phase, the CNN is used after the ultrasound imaging procedure has been performed and the signal data recorded, but before the imaging software has reconstructed an ultrasound image. The CNN is used as part of 702 of flowchart 700, finding the positions of the transducers, which is part of 312 of flowchart 300, reconstructing the image. The imaging procedure is performed at 308 of flowchart 300, with the signals recorded at 310 of flowchart 300 during the imaging procedure. The images are generated at 412-420 of flowchart 400, which is also part of 312 of flowchart 300, but done only after flowchart 700. Reconstructing the image, 312 of flowchart 300, is normally not done in real time immediately after the procedure, 308 and 310 of flowchart 300, but only some time afterwards, because the software may take a long time to reconstruct an image. In the application phase, the CNN produces as output precise positions of the transducers at the times the corresponding ultrasound signals were recorded from the subject. These precise positions of the transducers can then be used by the imaging software, starting with the rest of flowchart 70θ after 702, to generate accurate images and other data about the field of view of the subject being examined. In the application phase, for a given imaging procedure, the CNN may be used independently for each of a large number of partly overlapping subsets of the transducers, which optionally together cover the entire array of transducers. Each time the CNN is applied to input signals recorded from one subset of transducers, it produces output with the relative positions of the transducers in that subset. The output produced by the CNN using data from for each of the subsets of transducers produces information on the relative positions of transducers in that subset, optionally as a function of time if the input signals used for that application of the CNN were recorded over a long enough period of time that the positions of the transducers would be expected to change significantly, for example by more than 2% or 5% or 10% or 20% of an ultrasound wavelength. That information, for all the subsets, can then be integrated to find the relative positions of all the transducers in the array, or at least of a widely and densely distributed set of the transducers that can be used to locate the other transducers by interpolation. For example, the relative positions of the transducers in each of the overlapping subsets of transducers is combined using a least squares method to find the relative positions of all the transducers in the subsets.

Even if each application of the CNN only finds the positions of transducers in that subset at a fixed time, because all of the input signals it uses were acquired over a very short time period during which the transducers would not be expected to move significantly, the signal acquisition time period for all the applications of the CNN may be much longer, and the positions of the transducers may change significantly during that longer time period. In that case, the CNN output using data from each subset of transducers is optionally associated with its own time stamp, and when the different outputs are integrated to find the relative positions of all the transducers or a densely distributed set of the transducers in the array, the time dependence of the different outputs is taken into account, and the relative positions of the transducers is found as a function of time. The time dependence of the relative positions of the transducers may also be taken into account if one or more applications of the CNN do find the positions of the transducers in that subset as a function of time. In both these situations, the software that integrates the output from the different applications of the CNN takes into account that discrepancies in the position of a given transducer, that is part of two different overlapping subsets, may be due not only to random errors, which can be reduced by using a least squares fit, but may also be due to the different applications of the CNN using signals acquired at different times, between which the transducer has moved significantly.

Optionally, each subset of transducers, or at least some of the subsets of transducers, consist of all the transducers located within a bounded area of the surface of the body part being imaged, with the area being only a small fraction of the total area of the surface of the body part, for example less than 25%, or less than 10%, or less than 5%. Those subsets will sometimes be referred to herein as “patches.” Additionally or alternatively, at least some of the subsets include only some of the transducers within a small bounded area of the surface. Additionally or alternatively, at least some of the subsets include transducers that are distributed over the whole surface or a large part of it, for example more than 50% of the area of the surface, rather than being limited to a bounded area that is only a small fraction of the area of the surface. In that case, since each subset generally has only a small fraction of all the transducers in the array, the subset optionally has only a small fraction of the transducers that are within the area that the subset extends over.

Figures 8A and 8B schematically show two overlapping patches of transducers on the same array that is wrapped around a body part, such as a trunk. Such patches of transducers are optionally used as subsets of the transducers whose signal data is used by the CNN to find the positions of the transducers, as described above. They are referred to here as patches because they each consist of all the transducers in a small bounded area of the surface of the body part. What is described for the patches shown in FIGS. 8 A and 8B is also generally true for other subsets of transducers used for the CNN, even if they do not consist of all the transducers in a small bounded area of the surface of the body part. It should be understood that in a real system, the transducers will normally be much smaller and closer together relative to the dimensions of the array and the body part than in FIGS. 8 A and 8B, and there will be many more transducers in the array, and there will be many more transducers in each patch and/or the patches will be much smaller relative to the dimensions of the array and the body part.

In some embodiments of the invention, instead of the CNN using the signal from every combination of source-receiver pairs from a set of source transducers and receiver transducers that is in a patch, a more complicated set of source-receiver pairs may define the signals used by the CNN. For example, the CNN may use some receivers, but not other receivers, for signals transmitted by one source, and may use a different set of receivers for the signals transmitted by a difference source. But in the examples given here for purposes of describing the operation of the CNN, we only consider the case where the CNN uses signals from every combination of source transducer and receiver transducer, out of a set of source transducers and receiver transducers that is defined as belonging to one patch.

Figure 8A shows a patch 800 of a transducer array, with source transducers 802, 804, and 806, symbolized by triangles, and a larger number of receiver transducers 808, symbolized by shaded circles, making up the patch. The other transducers in the array, not part of this patch, are symbolized as unshaded circles. Figure 8B shows another patch 810 on the same array, seen from the same point of view, that partly overlaps patch 800. For example, source transducer 804 in patch 800 is the same as source transducer 812 in patch 810, and source transducer 806 in patch 800 is the same as source transducer 816 in patch 810. But source transducer 802 in patch 800 is not part of patch 810, and source transducer 814 in patch 810 is not part of patch 800. Similarly, receiver transducers 818 in patch 810 include some of the receiver transducers 808 in patch 800, but patch 81θ also has receiver transducers 818 that are not part of patch 800, and patch 80θ also has receiver transducers 808 that are not part of patch 810. In addition to patches 80θ and 810, the array has enough other patches, each of them overlapping somewhat with neighboring patches, to cover the whole array.

Each patch is associated with a set of recorded ultrasound signals, each of the signals transmitted from one of the source transducers in the patch to one of the receiver transducers in the patch. The set of recorded signals associated with a patch need not include all possible combinations of source transducers and receiver transducers that belong to the patch, but it may include all the possible combinations. Optionally, for at least some patches, the set of signal data associated with that patch is recorded over a relatively short time interval, compared to the total time it takes to record the signal data for the whole array. The sets of signal data associated with different patches are optionally recorded at different short time intervals, optionally not overlapping. This has the potential advantage that the positions of the transducers in each patch, to be found by the CNN, may not change very much during the time when the input signals to the CNN were recorded for that patch, so the CNN may not need to find the positions of the transducers as a function of time, which may require less computation. And even if the CNN does find the positions of the transducers as a function of time, the positions may be more accurate if the transducers do not move very much. For example, if there are 60 patches, and it takes 1 minute to record signal data for all the patches, and the signals, associated with a given patch, are recorded at various times that are distributed over the entire 1 minute, then the positions of the transducers might be expected to move significantly while the data is being recorded for that patch. But if all the signal data for any one of the 60 patches is recorded within 1 second, and the data for different patches is recorded consecutively, then we might expect the positions of the transducers not to move very much during the recording of signal data for each patch. In practice there is likely to be some overlap between the times of signal recording for different patches, if there are some sources that are included in more than one patch, since all of the signals transmitted by a given source are usually recorded at the same time. But as long as all the data for a given patch is recorded in a relatively short time, for example not much more than one second, then the CNN may give reasonably accurate results for the positions of the transducers for that patch.

Optionally, the absolute time of transmission from each source transducer, for example the clock time, is recorded for each recorded signal, and is included in the input to the CNN. Optionally, when the transducer position data found by the CNN for each patch, with or without time dependence, is integrated to find the relative position data of all the transducers, the software has access to the different transmission times for the different signals, and uses that information to estimate how the positions of the transducers were changing over time, during the time that all the signal data was recorded. Such changes in position of the transducers are expected to occur if recording all the signal data takes more than a few seconds, for example due to breathing, or even for a shorter time due to the subject’s cardiac cycle, or due to motion by the subject. Optionally, the software that integrates the output from the different applications of the CNN, using data from the different patches or subsets of transducers, uses models of the breathing cycle, the cardiac cycle, and/or likely motions by the subject, and how they affect the positions of different transducers, in order to better determine the positions of the transducers as a function of time.

Alternatively, especially if the positions of the transducers are not expected to change significantly during the entire time when signals are recorded, signals used as input to the CNN for a given patch may be recorded at times that are widely distributed over the entire time when signals are recorded.

When recording ultrasound signals, optionally only one source transducer transmits a signal at a given time, so that it is possible to separately record the signal for each source-receiver pair. In general signals received by many different receiver transducers from the same source transducer may be recorded separately at the same time.

In an exemplary transducer array that wraps around most of the surface of the trunk and is about 150 cm in circumference around the trunk and 40 cm in length along the trunk, there are, for example, about 24,000 transducers in the array, spaced about 5 mm apart, including, for example, about 6000 transducers that can act as source transducers, spaced about 1 cm apart. The spacing between transducers depends on the shape and size of the subject’s body, if the array is stretchable and flexible, and it should be understood that the spacing distances given here are only approximate examples. The width and length of the array may also vary depending on the shape and size of the subject’s body, and there may be different numbers of transducers used for different subjects.

This array is optionally designed to produce an image with 2 mm resolution, using ultrasound at 1 MHz with wavelength about 1.5 mm in soft tissue. Only one of the source transducers transmits ultrasound at a given time, and at that time the other transducers optionally all serve as receiver transducers. Alternatively, for each source transducer, only every other transducer is used as a receiver transducer, so there are 12,000 receiver transducers for each source transducer. There would then be a total of 72 million source-receiver pairs, each one producing a recorded signal that is used as input for generating an image. In an exemplary embodiment of the invention, each patch whose data is used by one application of the CNN comprises an area 10 cm by 10 cm, with 400 receiver transducers and 100 source transducers, and the patches are arranged on the frame of transducers as a Cartesian array. If there were no overlap between the patches, then there would be a total of 60 patches over the frame, 4 patches in the longitudinal direction by 15 patches in the azimuthal direction. If adjacent patches overlapped by 5 cm around the trunk and by 5 cm along the trunk, then there would be a total of 210 patches, 7 patches in the longitudinal direction by 30 patches in the azimuthal direction, and the CNN would be applied 210 times for each image, once using data from each patch. The patches need not be arranged in a Cartesian array, but can be arranged in any regular or irregular configuration. Optionally, some or all of the patches have dimensions that are less than 5% of the largest diameter of the field of view, or between 5% and 10% of the largest diameter, or between 10% and 20% of the largest diameter, or between 20% and 50% of the largest diameter, or more than 50% of the largest diameter. Optionally some or all of the patches have less than 20 transducers each, or between 2θ and 50 transducers, or between 5θ and 100 transducers, or between 10θ and 200 transducers, or between 20θ and 500 transducers, or between 50θ and 1000 transducers, or between 100θ and 2000 transducers, or between 200θ and 5000 transducers, or between 500θ and 10000 transducers, or more than 10000 transducers. Optionally, different patches have dimensions and/or numbers of transducers that are in different ones of these ranges.

Figure 9A shows a diagram 900 of the architecture of the CNN, according to an exemplary embodiment of the invention. Although other architectures can be used for the CNN, there is a potential advantage to using the U-shaped architecture of diagram 900, an architecture described by Ronneberger et al, “U-Net: Convolutional Networks for Biomedical Image Segmentation,” arXiv: 1505.04597vl [cs.CV] 18 May 2015. In this architecture, in the first part of the CNN, shown on the left side of FIG. 9A, the data set contracts at successive layers, using pooling operations, while the number of channels increases. The architecture of the contracting part of the CNN is similar to the architecture used in CNNs for image recognition, that determine, for example, whether the image as whole depicts a certain object, for example a cat or a dog. Successive layers integrate information from the image at larger and larger scales. But Ronneberger et al found that for segmenting a biomedical image, which requires assigning a class label to each pixel, it is useful to also have an expanding second part of the CNN, shown on the right side of FIG. 9A. In the layers in the expanding part of the CNN, upsampling is used to increase the size of the data set in successive layers, while decreasing the number of channels by merging channels. The right side of the CNN is more or less a mirror image of the left side, giving the CNN a U-shaped architecture. For different layers on the right side, data from a corresponding layer on the left side, with a data set of the same dimensions and the same number of channels, is copied and combined with the data on the right side, as shown by the three horizontal dashed lines in FIG. 9A. Ronneberger found that this high resolution information from the left side allows the layers on the right side to assemble a more precise output, identifying local structures in the image, much more efficiently than other CNN architectures that have been used for biomedical image segmentation. The inventor had the insight that the same U-shaped architecture would potentially be useful in a CNN for determining the positions of transducers from signal data, because the problem that must be solved has similarities to image segmentation. The recorded signal data from pairs of source and receiver transducers at known approximate positions is analogous to the image being segmented in Ronneberger et al. Assigning a precise position to each transducer is analogous to assigning a class label to each pixel in the image in Ronneberger et al. In both cases, local data is combined with more global data, to solve the problem.

A CNN algorithm with U-shaped architecture is also sometimes referred to as an encoderdecoder algorithm, as described for example in a web page of the Gabor Meili Research Knowledge Base, downloaded from <www(dot)gabormelli(dot)com(slash)Convolutional_(CNN/CNN)- based_Encoder-Decoder_Neural_Network>, on October 28, 2023. Ronneberger et al, cited above, refers to these algorithms as U-nets. The terms “U-shaped CNN,” “CNN with U-shaped architecture,” “U-net,” and “Encoder-decoder neural network” are used interchangeably herein.

In FIG. 9A, element 902 is a primary input data volume, sometimes referred to herein as the main input data volume, element 930 is a secondary input data volume, and element 940 is an output data volume of the whole CNN. The other numbered elements, in FIG. 9B as well as in FIG. 9A, all represent layers of the CNN, or in the case of element 932 a group of such layers. A layer of a CNN is a process that transforms a data volume into another data volume. In general, in a CNN, there is a series of layers, interconnected in a certain way, possibly branching or joining together, which perform a series of processes, starting with the one or more input data volumes, and at the end generates the output data volume of the whole CNN. In the description of FIGS. 9A and 9B, it should be understood that we sometimes informally use the label of a layer to refer to the data volume that results from that layer. For example, “data volume 928” means the data volume that results from layer 928, or is the output of layer 928, and similarly for the other layers of the CNN. Sometimes in the description of FIGS. 9A and 9B, a layer is not labeled explicitly or referred to as a layer at all, but its action is described, and it is only labeled implicitly by the label of the data volume that results from the layer.

If a U-shaped CNN is going to be used, because of the potential advantages described above, rather than a different type of CNN or machine learning algorithm, there are many different ways to construct a U-shaped CNN. Diagram 900 is one example of such an architecture, that is being used in tests by the inventor. Although the details described herein are not essential for the CNN to work well, certain overall features of the CNN described in diagram 90θ are believed to be potentially advantageous. The U-shaped structure is potentially advantageous, for the reasons given above. In addition, module 932 of diagram 900, which combines secondary input data in data volume 930 with the output 928 of the U-shaped part of the CNN operating on the main input data volume 902, is potentially advantageous, if secondary input data are used, in this case initial estimates of the positions of the transducers, obtained for example by modeling how the frame of transducers sits on the subject’s body, and optionally also from strain and/or angle sensors in the frame of transducers. Alternatively, such secondary input is not used, or it is combined with the data volume 928 in a different way than shown for module 932. The term “data volume” will be used herein to refer to a data array with any number of dimensions, not necessarily with three dimensions. In describing the CNN, the terms “data volume” and “data array” may be used for an input data volume, an output data volume, and any intermediate data volume resulting from processing by any of the layers in the CNN.

As used herein, an algorithm having a U-shaped architecture means an algorithm, such as the CNN algorithm outlined in FIG. 9A, that includes a contracting part, for example 906, 908, 91θ and 912 in FIG. 9A, followed by an expanding part, for example 914, 916, 92θ and 924 in FIG. 9 A. In the contracting part of the algorithm, an input data volume contracts, for example using pooling operations, in successive layers of the algorithm, while the number of channels increases. In the expanding part of the algorithm, the size of the data volume increases in successive layers, for example using upsampling, while decreasing the number of channels by merging channels. For at least some layers in the contracting part, a data volume is copied and combined with a data volume of the same dimensions and same number of channels in a corresponding layer in the expanding part. As used herein, a data volume is said to be passed through the U-shaped architecture of an algorithm when the data volume undergoes this process of successive contractions in the contracting part of the algorithm followed successive expansions in the expanding part of the algorithm, with data volumes in one or more layers of the contracting part being copied and combined with data volumes in corresponding layers of the expanding part.

Data volume 902 of diagram 900 is the main input data, which consists of recorded signals transmitted from each of a set of Ns source transducers in a patch, to each of a set of NR receiver transducers in the patch, each recorded signal having Nt values, one at each of a discrete set of times b. For example, the signal is recorded with a sampling time interval 6t, so ti = (i- 1 )5t. The sampling time interval, as noted above, is much shorter than an ultrasound wave period, for example by at least a factor of 2, or at least a factor of 5, or at least a factor of 10, or at least a factor of 20. Each piece of signal data is us,R(ti) for a source S in the patch, a receiver R in the patch, and a time b in the discrete set of times. Optionally, there is no need to list identifiers of the source and receiver transducers as part of the signal data in data volume 902, because, for example, the signal data may always be arranged in order of source transducers, and for each source transducer in order of receiver transducers, according to the order of source and receiver transducers, with their identifiers, listed in the secondary input data in data volume 930. Alternatively, identifiers of the source and receiver transducers are included as part of the signal data in data volume 902. The times in the signal data are measured, as noted above, relative to a defined time of transmission of the signal by the source transducer, for example, relative to the starting time of transmission by the source transducer, or relative to halfway between the starting time of transmission and ending time of transmission by the source transducer. The defined time of transmission is, for example, defined in the same way for all signals.

Optionally, in addition to using the raw signal data us,R(t) as input for the CNN, one or more attributes of the signal data are calculated from the signal data, and are included, together with the raw signal data, as part of input data volume 902. The number of attributes in input volume 902 is referred to as N_a, and there are N_a numbers in input volume 902 for each input signal US,R at each sampled time t. Using this additional data may improve the performance of the CNN. Examples of such attributes that can be calculated from the signal data and included in input volume 902 include the complex trace of the signal, defined in the paper by Taner et al cited above, the envelope of the signal, defined for example as the absolute value of the complex trace, and the instantaneous phase, defined for example as the phase of the complex trace. Other possible attributes, which are used sometimes in seismic imaging, include the signal with a frequency filter applied to it, time derivatives of the signal, and time integrals of the signal.

In general, input data volume 902 is a four-dimensional array of dimensions Ns x NR x Nt x N_a, and the data undergoes a convolution operation, using, for example, a 5 x 5 x 5 x N_a kernel, where the 125N_a values of the elements of the kernel are weights that are learned during the training of the CNN, as will be described in FIGS. 10A and 10B. The dot product of the kernel with successive 5 x 5 x 5 x N_a sub-arrays of the input data becomes part of the output data volume produced by layer 904 of the CNN. The data volume that results from layer 904 is informally referred to as “data volume 904.” The stride for choosing the different 5 x 5 x 5 x N_a sub-arrays is 1, for example, so that all the 5 x 5 x 5 x N_a sub-arrays are used, and the output of layer 904 would then be a data volume of dimensions (Ns-4) x (NR-4) X (Nt-4) x 1. Alternatively, before performing the convolution operation, the input volume is padded, for example with zeroes, for example with 4 extra elements in each of the Ns, NR, and Nt dimensions, two elements on each side for each dimension, and in that case the output would be a volume with dimensions Ns x NR x Nt x 1. Note that, because the kernel has the same size N_a in the attribute dimension as the input volume, and the input volume is not padded in that dimension, the output data volume is only one unit wide in the attribute dimension, so it is a 3-D array, not a 4-D array like the input volume. Alternatively, a kernel of different dimensions is used, and/or a different stride is used, and optionally a different amount of padding is used to keep the output data volume of the same dimensions as the input data volume.

This procedure is performed multiple times, for example 16 times, in general using a different kernel each time. Each time the convolution is performed, the output goes into a different channel of output data volume 904, which, in the case of 16 channels, and assuming the input volume was padded in the Ns, NR and Nt dimensions, has dimensions 16 x Ns x NR x Nt. Optionally, a bias, also a weight that is set by the training, is added to each element of the output data volume, optionally a different bias for each channel. Optionally, the elements of the output volume are activated, for example using rectified linear unit (ReLU) activation, in which negative values are replaced by zero. Alternatively, another activation method is used, such as tanh activation, replacing each element x by tanh(x), or sigmoid activation, replacing each element x by 1/(1 + e^-x), but ReLU activation has the potential advantage that it requires less computation. Optionally, the elements of the output data volume are normalized using one of several known methods of normalizing outputs in neural networks, for example batch norm, weight norm, layer norm, instance norm, group norm, batch-instance norm, and switchable norm. “Normalization techniques in deep neural networks,” by Aakash Bindel, Feb. 10, 2019, available at <medium(dot)com/techspace-usict/normalization-techniques-in-deep-neural-networks- 9121bfl00d8>, describes these normalization methods and variations on them, some advantages and disadvantages of each method, and some advantages of using normalization in neural networks.

In diagram 900, all of the data volumes that are outputs of a convolution operation, or outputs of a “dense” operation as will be described below, are activated with ReLU, and are normalized using instance normalization unless a different normalization method is specified below. But optionally, the activation and/or normalization methods, like other details of the CNN, are different, and are not necessarily critical to the performance of the CNN. In the next layer of the CNN, layer 906, data volume 904 then undergoes a convolution operation using a kernel, for example a l6 x 5 x 5 x 5 kernel, to produce a single channel of an output data volume, with all of the options described above for primary input data bloc 902. This is optionally repeated, with a generally different 16 x 5 x 5 x 5 kernel, for each channel of the output data volume, for example for 16 channels. Assuming that the input, data volume 904, is padded to keep the dimensions of the arrays the same, the output, data volume 906, again has dimensions 16 x Ns x NR x Nt.

A maximum pooling operation is then done, reducing the time dimensions by a factor of 2, replacing sets of two consecutive data elements in the time direction by the maximum of the two elements. Alternatively, a different pooling operation is done, for example taking an average of the two elements instead of the maximum. The number of channels is optionally doubled, for example from 16 channels to 32 channels, duplicating the data in each channel. Data volume 906 now has dimensions 32 x Ns x NR x (Nt/2). Data volume 906 then, in layer 908, undergoes a convolution operation, using a kernel, for example a 32 x 5 x 5 x 5 kernel, with any of the options described above, including using padding, adding a bias term, and using activation and normalization, to produce each channel of an output data volume 908, and this is repeated for each channel of data volume 908, for example 32 times with 32 generally different 32 x 5 x 5 x 5 kernels, to produce 32 channels of data volume 908, which will have dimensions 32 x Ns x NR x (Nt/2).

A pooling operation, for example a maximum pooling operation, is then done, in layer 910, to data volume 908, similar to the maximum pooling operation described for data volume 906, reducing the dimensions in the time direction to Nt/4, and the number of channels is optionally doubled again, to 64, duplicating the data in each channel. Data volume 908 then undergoes a convolution operation, using a kernel, for example a 64 x 5 x 5 x 5 kernel, to produce each channel of the output data volume, for example 64 generally different kernels to produce 64 channels of the output data volume, optionally adding bias terms and using activation and normalization as described above. The output is data volume 910, with dimensions 64 x Ns x NR x (Nt/4) for each attribute.

In layer 912, a pooling operation and a convolution operation similar to that described for data volume 908, with any of the options described there, are performed on data volume 910, increasing the number of channels to 128, and decreasing the time dimensions to Nt/8. The output is data volume 912, with dimensions 128 x Ns x NR x (Nt/8), assuming padding is used in the Ns, NR and Nt dimensions.

A convolution operation is performed, in layer 914, on data volume 912, with any of the options described above for convolution operations, but without a pooling operation and without changing the number of channels. The output is data volume 914, which, if padding is used in the Ns, NR and Nt dimensions, has the same dimensions at data volume 912.

Data volume 914 is then upsampled, doubling the time dimension to Nt/4, for example by duplicating each data element. The upsampled data volume 914, with dimensions 128 x Ns x NR x (Nt/4), is then concatenated with data volume 910, which has dimensions 64 x Ns x NR x (Nt/4), as shown by arrow 918, producing a data volume with 192 channels, with dimensions 192 x Ns x NR x (Nt/4). The number of channels is then decreased to 64, by performing a convolution operation on the data volume with 64 generally different kernels, one kernel for each of the 64 output channels, each kernel having dimensions 192 x 5 x 5 x 5. A convolution operation is then applied to this 64 x Ns x NR x (Nt/4) data volume, with any of the options described above for convolution operations. The output is data volume 916.

Data volume 916 is then upsampled, doubling the time dimension to Nt/2, for example by duplicating each data element. The upsampled data volume 916, with dimensions 64 x Ns x NR x (Nt/2), is then concatenated with data volume 908, which has dimensions 32 x Ns x NR x (Nt/2), as shown by arrow 922, producing a data volume with 96 channels, with dimensions 96 x Ns x NR x (Nt/2). The number of channels is then decreased to 32, by performing a convolution operation with 32 generally different kernels, one kernel for each of the 32 output channels, each kernel having dimensions 96 x 5 x 5 x 5. A convolution operation is then applied to this 32 x Ns x NR x (Nt/2) data volume, with any of the options described above for convolution operations. The output is data volume 920.

Data volume 920 is then upsampled, doubling the time dimension to Nt, for example by duplicating each data element. The upsampled data volume 920, with dimensions 32 x Ns x NR x Nt, is then concatenated with data volume 906, which has dimensions 16 x Ns x NR x Nt, as shown by arrow 926, producing a data volume with 48 channels, with dimensions 48 x Ns x NR x Nt. The number of channels is then decreased to 16, by performing a convolution operation with 16 generally different kernels, one kernel for each of the 16 output channels, each kernel having dimensions 48 x 5 x 5 x 5.A convolution operation is then applied to this 16 x Ns x NR x Nt data volume, with any of the options described above for convolution operations. The output is data volume 924.

Data volume 924 undergoes a convolution operation, with any of the options described above for convolution operations, without changing the Ns, NR and Nt dimensions of the data or the number of channels. The output is data volume 928.

Data volume 930 is a secondary input data volume, which contains initial estimates of the positions of the transducers in the patch. To find the initial estimates of position of each transducer, the following exemplary method is optionally used. The frame coordinates of each transducer are obtained, for example from the design specification of the frame of transducers, as described above in the description of FIG. 2. A mechanical model of the subject’s body is chosen, for example a model that is chosen to match the subject’s height, weight, body build, chest, waist and hip measurements, age, gender, and/or other characteristics. The body models available may be based on the population that is being served by the ultrasound system. This mechanical body model, whose purpose to calculate where on the body a given transducer will be located when the array is wrapped around the body, will generally be different from the body model described below in the description of FIG. 10A, whose purpose is to simulate ultrasound propagation in the body in order to train the CNN. The mechanical body model optionally includes the coordinate system used to reconstruct the image. Using known elastic properties of the frame, for example the elastic model of the array that may be included in the design specifications of the array, the frame coordinates of the transducers, the body model, an assumption about the posture of the subject during the data acquisition, for example lying prone on an examining table, and a model for how the frame is wrapped around the body part of the subject, initial estimates are calculated for the positions of the transducers in the coordinate system used to reconstruct the image.

Optionally, data from strain sensors and/or angle sensors incorporated in the frame is used to make the initial estimates of the transducer positions more accurate, by measuring distances and/or angles on the frame instead of relying on the body model and the elastic model of the array to calculate them. Optionally, the data from strain sensors and/or angle sensors is taken at the time of transmission of the ultrasound signal from each source transducer in the patch, and that data is used to calculate more accurate initial estimates of the transducer positions, for each of the signal transmission times. These initial estimates optionally include initial estimates of the position of each receiver transducer at the time of transmission of each source transducer, and an initial estimate of the position of each source transducer only at the time of its own signal transmission. Each of these initial estimated positions corresponds to one of the precise positions that the CNN will find as its final output. Alternatively, the strain sensor and/or angle sensor data is taken at fixed times throughout the period signal transmission for this patch, and recorded, and the recorded data, interpolating in time if necessary, is used to make more accurate initial guesses of the positions of the sensors at the time of signal transmission of each source transducer. Alternatively, for example if the transducers are not expected to change position significantly during the duration of signal transmissions for this patch, the strain sensor and/or angle sensor data is taken only once during the signal transmission for this patch, and is used to make a single improved estimate of the position of each transducer, which may be used in secondary data volume 930 for the initial estimates of the transducer positions for all times of transmission of the source transducers for this patch.

For each transducer, there are optionally 4 data elements, corresponding for example to the approximate r, θ, and z coordinates of the transducer, and an identifier for the transducer. In some embodiments of the invention, there is also a fifth data element, a time stamp, for each transducer, or even more data elements, as will be described below. But the main description of the CNN will describe the case where there are only 4 data elements for each transducer, and it should be understood that the description can be modified in a straightforward way if more than 4 elements are used. It is potentially advantageous for identifiers for the transducers to be included in the output of the CNN, because afterwards the CNN may be run using data from each of many different partly overlapping patches of transducers, and the output from the different applications of the CNN, using data from different patches, may be combined to find precise positions for all the transducers in the transducer array. The identifiers of the transducers will provide information on which transducers in the different patches are the same transducer, when the outputs for the different patches are combined. Even if one or more of the source transducers is also used as a receiver transducer when a different source transmitter is being used, optionally in that case the data for that transducer is entered twice in the secondary input data volume, once as a source and once as a receiver, so it will be known to the CNN that that transducer is used both as a source and as a receiver at different times.

Optionally, the recorded signal data in the main input data volume is entered in an order according to the order of the transducers listed in the secondary input data bloc. For example, if there were 3 source transducers and 16 receiver transducers in each patch, as shown in the simplified example of FIGS. 8 A and 8B, then in the secondary input data volume, the transducers are listed in the order SI, S2, S3; Rl, R2,...R16, and in the main input volume, the signals are listed in the order (S1,R1), (S1,R2),...(S1,R16), (S2, Rl), (S2, R2),...(S2,R16), (S3,R1), (S3,R2),...(S3,R16). Then the identity of the source transducer and the receiver transducer will be known for each recorded signal, and it will not be necessary to include that information with the recorded signal data in the main input volume. Alternatively, identifiers for the source and receiver transducer for each recorded signal are included with the recorded signal data.

Optionally, each of the Ns source transducers transmits a data signal simultaneously to all NR receiver transducers. Because the travel time of ultrasound across the body is very short, less than a millisecond, which is not enough time for the transducers to move significantly, the acquisition time of all the signals transmitted from a given source transducer is effectively the same for all the receiver transducers. The position of each source transducer is based on the signals transmitted from that source transducer to each of the receiver transducers, so there is only one position that the CNN will output for the position of each source transducer. But each receiver transducer receives signals from many different source transducers, received at different times. In the case where the CNN is taking time dependence into account, it will output Ns different positions for each receiver transducer, each position based on the signal it receives from a different one of the Ns source transducers. When the image is reconstructed using pre-stack depth Kirchhoff migration, the travel time of the signal from each source transducer to that receiver transducer will be based on the position of that receiver transducer at the time the signal was sent by that source transducer. Assuming the CNN produces accurate values for the transducer positions, the reconstructed image will be accurate, even if the transducers are moving by much more than a small fraction of an ultrasound wavelength during the data acquisition time. Then the CNN will produce 4NS(NR + 1) pieces of data as output, 4Ns pieces of data for the single positions of the Ns source transducers, and 4NSNR pieces of data for the Ns different positions (at different times) of each of the NR receiver transducers. So the output of the CNN should be an array with 4NS(NR + 1) elements. Optionally, secondary input data array 930 has the same dimensions as the output array, because of the architecture of module 932, as will be explained, and also because there might be a different initial estimate for the position of each receiver transducer at the time it receives a signal from each source transducer. Even though the body model will generally give only a single initial estimate for the position of each receiver transducer, data from the strain and/or angle sensors in the array could give more accurate initial estimates that are different at each time. If strain sensor and angle sensor data is not used, and there is only a single initial estimate for the position of each receiver transducer, then identical values can be used for the Ns initial estimates of position of each receiver transducer.

In some embodiments of the invention, the CNN never takes into account the changing positions of the receiver transducers over time. In such an embodiment, only one position is found for each receiver transducer, and the output of the CNN, and the secondary input data volume 930, may each have 4(Ns + NR) pieces of data, instead of 4NS(NR + 1) pieces of data. The dimensions of some of the intermediate data volumes in module 932, described below, are accordingly also changed. But in the description below, the CNN does take into account the changing positions of the receiver transducers, and the output and the secondary input data volume both have 4NS(NR + 1) pieces of data.

In some embodiments of the invention, in addition to using the recorded signal data and the initial estimates of the transducer positions, at the time of each transmitted signal, as input to the CNN, the time at which each signal is transmitted, for example the absolute clock time, is also used as input to the CNN, for example as part of secondary input data volume 930. Optionally, a recorded time stamp for each transmitted signal is used for this purpose, for example the time stamp that is described above as optionally recorded for the signal transmitted by each source transducer. For the purposes of the CNN the transmission time need only be known to within the time it takes for a transducer to change its position significantly, for example by 10% of an ultrasound wavelength. Optionally, for each receiver transducer for each signal it receives, the time of signal transmission is included as a fifth piece of data, in addition to the three coordinates (r, θ, z) of the initial guess of its position, and its identifier. Similarly, for each source transducer, the time of its signal transmission is included as a fifth piece of data, in addition to the three coordinates (r, θ, z) of the initial guess for its position, and its identifier. Then there will be 5NS(NR + 1) pieces of data in secondary input volume 930, instead of 4NS(NR + 1) pieces of data. Optionally, the time of transmission of each signal is included also in output data volume 940, at least because of the structure of the CNN, which has the same dimensions for secondary input volume 93θ and output data volume 940. In addition, including the time of transmission in the output of the CNN may be useful for the software that is integrating the outputs of the different applications of the CNN, corresponding to different patches of transducers, because the input to that software is the output from each of the applications of the CNN. That software may use information about the times of signal transmission in order to reconcile discrepancies in the output data for different patches, for example by modeling likely causes of change in transducer position, such as cardiac and breathing cycles.

Optionally, additional pieces of input data are used, for each source transducer, and for each receiver transducer at the time it receives a signal from each source transducer, resulting for example in 6NS(NR + 1) or 7NS(NR + 1) pieces of data in secondary input volume 930. In embodiments where the time of signal transmission, and/or other data, is included in secondary input volume 930 for each source transducer, and for each receiver transducer for the signal received from each source transducer, it should be understood that, in the description below, data volumes that are described as having dimension 4NS(NR + 1) would instead have larger dimensions, for example 5NS(NR + 1) or more, as the case may be. Alternatively, the time of signal transmission, and/or other additional data, are not included as an extra piece of data for each source transducer and for each receiver transducer at the time of signal transmission by each source transducer. For example, the time of signal transmission coming from a given source transducer is the same for each receiver transducer that receives it, so it is redundant to enter this same transmission time for every receiver transducer. In this case, signal transmission times may be included in secondary input data volume 930 for each source transducer, with the understanding that the signal transmission time is the same for all receiver transducers when they receive a signal from that source transducer. In that case, for example, the dimensions of secondary input data volume would be NS(4NR + 5), and it should be understood that in the description below, data volumes that are now described as having dimensions 4NS(NR + 1) could instead have dimensions NS(4NR + 5). Including the time of signal transmission in the input data explicitly, instead of only labeling it implicitly by the source transducers that transmitted it, has some potential advantages. For example, for finding the position of a receiver transducer at the time it receives a signal from a certain source transducer, it may be useful to know how much time has elapsed since it received signals from various other source transducers. That will affect how much the position of the receiver transducer is likely to have changed between one signal and another. So these times of signal transmission are potentially useful input data that could make the CNN more accurate in finding the positions of the transducers at the time of each signal.

Secondary input data volume 930 is combined with data volume 928, whih comes ultimately from primary input data volume 902, in a module 932. In some embodiments of the invention, data from the secondary input data volume, optionally after processing in one or more layers, is combined with data from the primary input data volume, optionally after processing in one or more layers, using one or more combining procedures, in parallel with each other, in series with each other, or both. Optionally, the two data volumes that are combined in one or more of the combining procedures have the same dimensions, and the two data volumes are combined element- by-element, producing an output data volume of the same dimensions as each of the two input data volumes of the combining procedure. For example the combining procedure comprises element- by-element addition, or element-by-element multiplication. Alternatively, the combining procedure comprises concatenation, and if the two input data volumes to the combining procedure have the same dimensions, then the output data volume has twice the dimensions of the input volumes. Optionally, there are at least three combining procedures, one of them being element- by-element addition, one of them being element-by-element multiplication, and one of them being concatenation. That is the case, for example, in the exemplary embodiment described herein for module 932.

The details of how secondary input volume 930 is combined with data volume 928 in an exemplary embodiment of the invention are shown in FIG. 9B, where a dashed rectangle surrounds the elements that make up module 932. Data volume 928 may be referred to herein as the first intermediate data volume, the result of passing primary input data volume 902 through the U- shaped architecture from 906 through 928. Data volume 928 undergoes a flattening operation, in which the 16 x NR x Ns x Nt four-dimensional array is replaced by a one-dimensional array with 16NRNSNI elements. That array then goes through a dense layer, in which every element of the output is a linear combination of every element of the input, to yield a one-dimensional data volume 931 with 4NS(NR + 1) elements. Meanwhile, secondary input data volume 930 goes through a dense layer 933, in which every element of the output is a linear combination of every element of the input, yielding a one-dimensional data volume 933 with 4NS(NR + 1) elements. Data volumes 931 and 933 are then added together element by element, producing another one-dimensional data volume 934 with 4NS(NR + 1) elements. Data volume 934 may be referred to herein as the second intermediate data volume, the result of combining data volumes 928 and 930 by element-by- element addition after first passing both of them at least through dense layers. Data volume 934 then passes through a dense layer 935, yielding a one-dimensional data volume 935 with 4NS(NR + 1) elements. Data volume 931 and data volume 935 are then then multiplied together, element by element, yielding a one-dimensional data volume 936 with 4NS(NR + 1) elements. Data volume 936 may be referred to herein as the third intermediate data volume, the result of combining data volumes 928 and 934 by element-by-element multiplication after first passing both of them at least through dense layers. Data volume 936 in then concatenated with data volume 930, to yield a onedimensional data volume 937 with 8NS(NR + l)elements. Data volume 937 may be referred to herein as the fourth intermediate data volume, the result of combining data volumes 936 and 930 by concatenation.

Returning to FIG. 9A, data volume 937, the output of module 932, goes through a dense layer 938, producing a data volume 938 with 16NS(NR + l)elements. Data volume 938 then goes through another dense layer 940, producing a final output data volume 940 with 4NS(NR + 1) elements. These elements are supposed to correspond to the precise r, θ, and z coordinates, and the identifier, for each the NR transducers at the time it receives a signal from each of the Ns source transducers, and the precise r, θ, and z coordinates and the identifier for each of the Ns source transducers in the patch.

The precise coordinates for each of the transducers in the frame are optionally found by running the CNN for a large number of overlapping patches, and running software that combines output data from different patches. Optionally, this software takes into account that the positions of transducers might have changed during the data acquisition time for different source transducers and different patches, even if the CNN does not take changing transducer positions into account. If the CNN outputs give different values for the position of a transducer that is part of two different overlapping patches, then the software that combines the outputs reconciles the difference by assuming that it might be due to changes over time, or it might be due to a random error.

Method of training the convolutional neural network

The convolutional neural network, such as the CNN described in FIG. 9A, is trained using a machine learning algorithm, before it is used to find precise positions of the transducers. The training, for example, finds values for the weights that are to be used for each channel each time a convolution operation is performed, as well as each time a dense layer is used, as well as any other weights that are used by the CNN. An algorithm for generating a very large number of training examples is defined, as described in FIG. 10A. A very large number of the training examples, for example training examples chosen at random from the full space of training examples, are used sequentially to train the CNN, as described in FIG. 10B. Choosing the training examples at random from the space of possible training examples, or otherwise ensuring that the training examples are well distributed around the full space of possible training examples, is potentially advantageous because it tends to ensure that the weights found will be weights that are optimized for the full space of parameters that the CNN is going to be used for, rather than being weights that are optimized for only a small part of the parameter space. Training is done first with a large number of training examples to generate a set of weights that works well, before the CNN is used for the first time in production runs to generate images, for example in element 702 of FIG. 7A. Optionally, training continues with further training examples, even after an effective set of weights is generated and production runs have been done, to further improve the weights over time.

Block diagram 1000 of FIG. 10A shows a method of generating training examples that are expected to sufficiently cover the space of transducer locations, body models, and patches that the CNN will encounter in the application phase, so that the CNN works well with almost all the input data from subjects that it encounters in the application phase. At 1002, several body models, or body part models, are generated, based for example on the population of subjects that the CNN is intended to be used for. Each body model or body part model includes ultrasound propagation properties, including for example sound speed, acoustic impedance, and attenuation, optionally for both p-waves and s-waves, at least at the frequency that is going to be used for the CNN, or a function of frequency for a range of frequencies that is going to be used, all as a function of location (r, θ, z) in the body model, and a body surface shape, giving the shape of the outer surface of the body or body part on which the transducers are located. Each body surface shape corresponds, for example, to a subject of a particular age, gender, and body type. Each model of ultrasound propagaition properties corresponds, for example, to a different distribution of bones and soft tissue within the body, as well as to different values and distributions of sound speed, acoustic impedance and attenuation, for p-waves and s-waves, among the bone voxels and among the soft tissue voxels. Optionally, there are relatively few different body models generated at 1002, for example 20 or 50 or 100 or 200 or 500 different body models, or a larger, smaller, or intermediate number, while finer gradations in body models are covered by the variations that will be generated from each body model. Optionally, the body surface shape and distribution of tissue types for some or all of the body models are based on medical images, for example MRI images, of people belonging to the population of potential subjects.

At 1004, small variations are generated for the ultrasound propagation properties model for each body model. For example, each body model has a different spatial distribution of tissue types, such as bone and soft tissue, and different types of bone and different types of soft tissue, with different sound speeds, acoustic impedance and attenuation, for p-waves and s-waves. In addition, each variation of the ultrasound propagation properties model might have different sound speeds, acoustic impedance, etc., for a given tissue type. The variations in boundaries of the different tissue types, and the variations in ultrasound propagation properties for the different tissue types, optionally reflect the distribution in 3-D shapes of different body parts, such as particular bones and organs, among the population of potential subjects, and the distribution of average ultrasound propagation properties in a given body part, among the population of potential subjects. It should be understood that it is not necessary, at 1004, to explicitly calculate a full body model for each variation. Rather, at 1004 an algorithm is established for generating such variations. For example, at each of a large number of boundary points between different tissue types, the algorithm specifies that the location of each boundary point can vary by ± 5 mm in units of 1 mm, and the variations in the locations of nearby boundary points are correlated in a realistic way. Similarly, for example, the algorithm specifies that the sound speed for a given tissue type varies by ±2.5% in units of 0.5% at different points inside that tissue type, with variation of sound speed as a function of spatial location at nearby points varying in a realistic way. For example, for one or both of the location of the boundary between tissue types, and the ultrasound propagation properties as a function of location within a given tissue type, cubic splines are used, to ensure that the boundary between tissue types is realistically smooth, and to ensure that the variation in ultrasound propagation properties as a function of position for a given tissue type is realistically smooth.

At 1006, for each body shape, a calculation is made of the expected body surface positions of transducers, when the flexible stretchable array of transducers is stretched over the body. Optionally, this calculation takes into account both the elastic properties of the array, and the deformability of the body surface by the forces generated by the array. This calculation may be similar to the calculation, described above in the description of the CNN in FIG. 9B, for generating initial estimates of the positions of the transducers on the subject’s body. At 1008, an algorithm for generating variations in the 3-D positions of the transducers is established, including variations in the r, θ, and z position of each transducer, due to variations in the exact shape of the body surface, and variations in the exact elastic properties of the array. Correlations between the displacements of nearby transducers are used in this algorithm. For example, the elastic modulus between each pair of adjacent transducers is allowed to vary by ±10%, in units of 2%. Optionally, for each training example generated at 1006 and 1008, one or more additional training examples are generated by squeezing, stretching, flipping or deforming the body model, including the positions of the transducers.

At 1010, an algorithm for generating patches is established. The patches are generated, for example randomly, in the same way that patches will be generated for the application phase of the CNN, for example patches of the same number of source transducers, the same number of receiver transducers, and covering the same size and shape of surface area, as the patches that will be used by the CNN in the application phase. Alternatively, the same set of patches is always used by the CNN, and only that set of patches is used in the training.

If the CNN is going to take into account the changing positions of receiver transducers at the different times of signal transmission for the different source transducers, then the training examples should also include changing positions of the transducers during the time period over which signals are transmitted, and the following additional step is used in 1008 to generate training examples for the CNN. For each patch, the set of source transducers is determined, and for each training example, a different time of signal transmission is specified for each source transducer. Optionally, the same signal transmission time is always specified for a given source transducer for a given patch, and these signal transmission times are always used when imaging a subject. Alternatively, the signal transmission times for each of the different source transducers are not always the same when imaging different subjects, but are selected according to some random or non-random procedure, and the same procedure is used to select signal transmission times for the training examples. For each of the signal transmission times, a different set of small variations is generated for the positions of all the transducers in the patch. Optionally, these small variations do not vary independently among transducers, and from one signal transmission time to the next, but are correlated in time and over the body of the subject in a realistic way, for example taking into account the breathing cycle, the cardiac cycle, and using a model for small movements that the subject might be expected to make, even when trying to remain still. If the transmission signal data is going to be gated to the cardiac cycle and/or the breathing cycle, then that is taken into account in modeling the effect of the cardiac cycle and/or the breathing cycle on the small variations in positions of the different transducers at the different signal transmission times. Alternatively, variations in transducer position over time and over the body of the subject are uncorrelated. However, including realistic correlations in the variations in the transducer positions in time and space might make the CNN work more accurately when it is trained. Box 1012 may be thought of as a collection of training examples that will be used at random for the training session that will performed in flowchart 1013 of FIG. 10B. However, it is not necessary to actually generate all the very large number training examples that could be used, before the training begins, which would be very wasteful of computer memory. Instead, box 1012 optionally has a set of rules for generating training examples, using the procedures described for 1002, 1004, 1006, 1008, and 1010. Each time a new training example is needed by training method 1013, it is generated by the rules of box 1012, and the details of that training example can be discarded from the computer memory before the next training example is generated.

At 1014 of FIG. 10B, a new training example is generated. At 1016, the 3-D wave equation is solved for each of the source transducers in the patch for this example, and for each source transducer the signal that would be received by each receiver transducer in the patch is calculated, thus generating a set of Nt data points for each source-receiver pair in this patch. If the CNN is being trained to take into account the changing positions of the transducers at different times during the transmission of signals by different source transducers, then the differing positions of transducers at the different times of transmission of signals, in this training example, are taken into account in calculating the signal that would be received by each receiver transducer from each source transducer. At 1018, the CNN is run, using the signals calculated in 1016 as the main input. For the secondary input to the CNN, positions of the transducers in this patch are generated, that resemble the initial estimated positions of the transducers that would be found from the body model and the elastic model of the array, optionally supplemented by data from strain and/or angle sensors in the flexible stretchable frame of the array. For example, starting with the exact positions for the transducers generated in 1008, error terms, generated by a pseudo-random number generator with a Gaussian distribution and optionally similar in magnitude to the precision of the strain and/or angle sensors, are added to the r, θ, and z coordinates of each transducer. Alternatively, a more precise simulation is made of the distribution of errors in position introduced by the sensors, for example taking into account the operation of the sensors and the algorithm used to calculate the positions of the transducers from the sensor data. Alternatively, any other method of introducing errors into the positions of the transducers is used, for example rounding off or smoothing the positions of the transducers generated in 1008.

Running the CNN, with the current set of weights, produces calculated 3-D positions for all transducers in the patch. If the CNN is taking into account the changing positions of the receiver transducers at the time of transmission of different source transducers, then output of the CNN will include a calculated 3-D position for each receiver transducer at the time of transmission of each source transducer in the patch, as well as a single position for the 3-D position of each source transducer at the time it transmits its signal. Meanwhile, at 1022, the exact positions of the transducers in this patch are looked up for this training example, as generated at 1008, including the different positions of the receiver transducers at different transmission times if the CNN is taking that into account. At 1024, the exact positions of the transducers in this patch, for this training example, are compared to the positions of the transducers calculated by output of the CNN at 1020, and a measure of the difference is found. It should be understood that initially, arbitrary starting values may be used for all weights in the CNN, for example all the weights are set at 1, and the actual positions of the tranducers at 1022 will have very poor agreement with the positions of the transducers according to the output of the CNN. At 1026, the measure of this difference is back propagated through the CNN. Because the layers of the CNN are generally more or less linear functions, back propagating the difference through the CNN provides information on how the weights should be changed, in order to reduce the difference between the actual positions of the transducers, and the positions calculated by the CNN. At 1028, the weights of the CNN are updated to values for which, according to this linear approximation and according to this training example, the actual positions of the transducers will be equal to the positions calculated by the CNN.

At 1030, an evaluation is made as the whether the agreement between the actual positions of the transducers in this patch, and the positions calculated by the CNN is good enough that it is not necessary to update the weights any further. Optionally, this is not done until a fairly large number of training examples have already been evaluated, and until the latest changes in the values of the weights, for example the relative changes, are very small for all the weights. If it is decided at 1030 that further adjustment of the weights is needed, then a new training example is chosen at 1014, for example at random using the algorithm for generating training examples in 1012. The training continues, typically for a very large number of training examples, that are well distributed over the space of training examples, until the weights converge to a set of values that gives good agreement between the actual positions of the transducers used in that training example, and the output of the CNN, and continues to give such good agreement for a large range of training examples. Once the weights have converged, they are outputted at 1032, and the fully trained CNN may be used for production runs. But as noted above, optionally further training takes place later, producing an even better set of weights.

How errors in velocity model affect Me._z(t, a)

Figure 11A schematically shows a cross-sectional view 1100 of a field of view of a body part that is being imaged by ultrasound. The same diagram could apply also to a field of view of a geological formation that is being imaged by seismic imaging. Figure 11 A illustrates why an error in the velocity model of the field of view results in a given feature in an angle gather appearing at a different depth beneath the surface for finite reflection angle a than for a = 0, and why the angle gather data can be used to correct an error in the velocity model.

In FIG. 11 A, the body surface 1102 is shown as flat, and an angle gather is being constructed for an image point 1104 beneath the surface of the body, with coordinates (r, θ, z). The r coordinate is a measure of depth beneath the body surface, and the angle gather begins at the surface of the body at a location 1106 directly above the image point, with the same coordinates (θ, z) as the image point. The velocity model is assumed to have a sound speed V that does not depend on the spatial coordinates r, θ, and z, so the ray optics paths of ultrasound waves are straight lines. The velocity model error that the angle gather data will reveal is an error in the average velocity between the surface of the body and the depth of the image point, and averaged laterally as well. It is not a measure of any lateral variations in the sound speed, which will be ignored here. The migrated reflectivity Me,_z(t, a), for a close to 0, is represented by the amplitude of an ultrasound wave that is transmitted from a source transducer 1108 very close to location 1106, reflecting off image point 1004, and received by a receiver 111θ also very close to location 1106 on the opposite side from transducer 1108. In the ray optics approximation, the sound wave follows straight line path 1112 from source 1108 to image point 1104, and a straight line path 1114 from image point 1104 to receiver 1110. The migrated reflectivity Me,_z(t, a) for a larger value of a, for example on the order of 45°, is represented by the amplitude of an ultrasound wave that is transmitted from a source transducer 1116 a larger distance x = r tana to the left of location 1106, and received by a receiver transducer 1118 located at the same distance x to the right of location 1106. The ultrasound wave from source 1116, in the ray optics approximation, follows a straight line path from source 1116 to image point 1104, and a straight line path 1122 from image point 1104 to receiver transducer 1118, reflecting from image point 1104 at a reflection angle a.

It should be understood that diagram 1100 is full of very rough approximations, but nevertheless it illustrates why an error in the velocity model V will result in an a-dependence of the migrated reflectivity Me,_z(t, a) that can be used to correct the error in the velocity model. For example, in general, the surface of a body part that is being imaged by ultrasound will not be flat, but will be curved, and will not be symmetric around the location 1106 of the top of the angle gather. The sound speed V will generally not be uniform as a function of r, θ, and z, and the ray optics paths of the ultrasound waves will be curved, rather than straight lines. The migrated reflectivity Me,_z(t, a) does not depend only on the signal from source 1108 to receiver 1110 for the case of a close to 0, and does not depend only on the signal from source 1116 to receiver 1118 for the case of tana = x/r, but in each case will be the sum of signals from many source-receiver pairs, and the receiver signal in each of those cases will come not just from that one source, but from many other sources as well, and the signals will not be just the result of reflecting from image point 1104, but from many other points in the field of view as well. However, the migrated reflectivity can be treated as if it is due only to signals from source-receiver pairs with similar geometry to source 1108 and receiver 1110 in the case of a close to 0, reflecting from image point 1104, and with similar geometry to source 1116 and receiver 1118 in the case of larger a, spaced at the same distance x = r tana on opposite sides of location 1106, reflecting from image point 1104. That is because signals from other source-receiver pairs, and reflecting from other image points 1104, may all tend to cancel out by destructive interference, which is the principle behind migration. So the migrated reflectivity Me,_z(t, a) for a given a can be treated as if it comes only from a source transducer and a receiver transducer located at a same distance x = r tana in opposite directions along the surface from location 1106.

The time needed for an ultrasound wave to travel from source 1108 to image point 1104, or from image point 1104 to receiver 1110, is to. This number does not depend on the velocity model being correct, since the echo time 2to from source 1108 to image point 1104 to receiver 1110 is measured directly. But if the velocity model is not correct, then the echo time will not be calculated correctly for ultrasound waves reflected from image point 1104, sent from a source to a receiver that are located off to the sides of location 1106, such as source 1116 and receiver 1118 that are each located a distance x to the sides. If the velocity V is correct, then image point 1104 is a distance r = Vto beneath the surface. In the general case, where V might vary with depth r, then it is an average velocity between the surface and image point 1104 that relates r to to, but in the simplified picture being considered here, V is uniform, and it is not necessary to distinguish between local velocity and average velocity. The ultrasound travel time from source 1116 to image point 1104, or from image point 1104 to receiver 1118, will be [to² + x²/V²]^1/2, assuming the velocity V is correct. Based on the assumed velocity V, the signal received at time [to² + x² /V²] ^1/2 by receiver 1118 from source 1116 will be assigned to Me,_z(to, a) where a = arctan(x/Vto). That is, the signal will be assigned to Me,_z(t, a) where t is the echo time times 1/(1 + x²/V²to²)^1/2. If the velocity model has a small error AV, with the true velocity being V + AV, then the true ultrasound travel time from source 1116 to image point 1104, or from image point 1104 to receiver 1118, will be [to² + x²/(V+AV)²]^1/2. The echo from source 1116 to image point 1104 to receiver 1118 will be assigned not to t = to, but to t equal to the true echo time times 1/(1 + x²/V²to²)^1/2. The error in the assigned value of t will be At = [to² + x²/(V+AV)²]^1/2/[l + x²/V²to²]^1/2 - to. To first order in AV, this is At = -to(AV/V)sin²a. Angle gathers, in which the reflectivity data is sorted by reflection angle a, are only one type of gather that can be used to correct the velocity model. For example, offset gathers, in which the reflectivity data is sorted by the offset distance X between the source and receiver transducer, can also be used to correct the velocity model. This can be done because, for a given depth r of the image point, or a given vertical echo time to, the offset distance X is generally related to the reflection angle a. An average velocity error AV will give rise to an error At that depends on the offset distance X, or on any other similar parameter that the gathers use to sort the reflectivity data. By estimating how At varies with a, X, or any similar gather parameter that characterizes how oblique ray paths differ from vertical ray paths, the velocity error AV can be estimated, as will be described below.

Figure 11B is a schematic illustration of two examples 1124 and 1132 of simplified migrated reflectivity data Me,_z(t, a), to show how an error in the velocity model affects the appearance of the data. Both examples are based on a simplified image that consists of only a single image point, such as image point 1104 in FIG. 11 A, that strongly reflects ultrasound sound, and nothing else in the field of view reflects ultrasound. The plots show vertical time t (the ultrasound wave travel time from the top of the angle gather, such as location 1106 in FIG. 11 A, down to the reflecting image point) on a vertical axis 1126, vs. reflection angle a on a horizontal axis 1128. For a given reflection angle a, the vertical travel time is calculated from the observed echo time following an oblique path, assuming a particular velocity V. In example 1124, the velocity V is assumed to be correct. Dots 1130, at each value of a, all show the same vertical time t, which is to, the vertical travel time for image point 1104 in FIG. 11 A. Because the velocity model is correct, the oblique travel times, at finite values of a, are all calculated correctly, and the vertical travel time, which only depends on the depth of the reflecting point beneath the surface, is the same for all values of a.

In example 1132, the velocity V is assumed to have an error of AV, with the real velocity being V - AV. As a result of this error in the velocity, the calculated vertical travel time has an error At = to(AV/V)sin²a, as calculated from FIG. 11 A. At a = 0, the vertical travel time is correctly given by to. But at higher values of a, the calculated vertical travel time is slightly greater, and the dots 1134 fall off in proportion to sin²a. This change in the appearance of the migrated reflectivity Mo,_z(t, a) can be used to determine errors in the velocity model at each angle gather (θ, z), and will be described in FIGS. 12 A and 12B.

As described in FIG. 7A, the correction in the velocity model, for a given angle gather (θ, z), is found by observing features in the migrated reflectivity Mo,_z(t, a), similar to the fall off in dots 1134 in FIG. 1 IB, using those features to estimate AV, correcting Me,_z(t, a) by replacing t by t + At, and then examining the corrected Me,_z(t, a) to see if the velocity model still has an error AV, and iterating until Me,_z(t, a) no longer implies a velocity error AV, for example if the dots in FIG. 11B go straight across like dots 1130, and the calculated AV vanishes. This iterative method of zeroing in on the velocity error means that it does not matter if the geometry shown in FIG. 11 A is very simplified and not realistic. As long as it gives an a-dependence of Me,_z(t, a) that has the right sign when there is a velocity error AV, we can find the velocity error.

Using the a-dependence of Me,_z(t, a) to find the velocity error AV, which is really an error in the average velocity V_rms(r, θ, z) defined above, ignores the dependence of velocity on θ and z, for example it does not take into account that oblique ray paths of the ultrasound, such as paths 112θ and 1122 in FIG. 11 A, traverse regions of the field of view where θ and/or z have values other than the values of θ and z for the angle gather being considered. To some extent, the effects of the θ and z dependence of V on the travel time of ultrasound waves can be ignored, because their effects are generally small compared to the effects of an error in the velocity at the θ and z of the angle gather. A uniform gradient of V with 0 or z will not have any effect on the total travel time from source to receiver in FIG. 11 A, because the error in the travel time from the source to the image point in FIG. 11A will be equal and opposite in sign to the error in the travel time from the image point to the receiver, and only second and higher derivatives of V with θ and z will matter. Those relatively small effects are taken into account in the outer iteration loop of FIG. 7A, where the velocity V(r, θ, z) is found from V_rms(r, θ, z) at element 736, smoothing over θ and z, thereby taking into account the θ and z dependence of the velocity, and the updated velocity V is then used to repeat the whole procedure, from 706 to 734, of evaluating Me,_z(t, a) for each angle gather (θ, z) and finding AV for each angle gather, iterating this outer loop until V(r, θ, z) converges globally.

Figure 1 IB shows how the velocity error AV could be evaluated in the artificial simplified case where there is only a single point in the field of view that reflects ultrasound. When generating a real ultrasound image, there are a large number of points of different reflectivities distributed throughout the field of view. Figures 12A and 12B illustrate how AV can be found from Me,_z(t, a) in such a case. Figure 12A shows a plot 1200 of M(t, a) for a particular angle gather. This data was actually taken from seismic imaging of a geological formation, but the method that it illustrates can be used also in ultrasound medical imaging. Time t, in milliseconds, is shown on vertical axis 1202, increasing going downward, and reflection angle a is shown on horizontal axis 1204. M is plotted as a function of t for each of 14 different values of a, ranging from 8° to 40.5° at intervals of 2.5°. For each of these values of a, the amplitude of M, in arbitrary units, is shown horizontally. To make it easier to compare levels of M at different times, the space under each M(t) curve is filled in black whenever M(t) is positive. Each M(t) curve oscillates irregularly with a varying period on the order of 25 milliseconds, which corresponds to the pulse length of the transmitted sound wave, and the spatial resolution of the image. M(t, a) has a high amplitude at those vertical travel times t for which the rock in the field of view has a high reflectivity of sound waves. The error in the velocity V is found by looking at M(t, a) at a specific time t, for example at t = 1110 ms, shown by horizontal line 1208. Figure 12B shows the amplitude M, on vertical axis 1214, as a function of reflection angle a, on horizontal axis 1212, at t = 1110 ms, taken from FIG. 12A. Small squares 1216 show the amplitude M at each of the values of a for which M(t, a) is plotted in FIG. 12A. These points have been fitted to the function of I + Gsin²a, and the best fit is shown as curve 1218. The fact that the best fit has non-zero G, resulting in a downward drift of M with increasing a, shows that there is a small error in the velocity model. We expect, from the analysis of FIG. 11 A, that the relative velocity error AV/V will be proportional to G/I. The actual expression, consistent with this, is derived by Swann, cited above, and is given above in the description of element 728 of FIG. 7A.

Pre-stack and Post-stack attributes

Figure 13 is a flowchart 1300 showing further details of how pre-stack and post-stack attributes are found in 41θ and 414 of FIG. 4. At 1302, the reflectivity M(r, θ, z, a), as a function of position (r, θ, z) and reflection angle a, is used to find a correct velocity model, giving the sound speed (interval velocity) V(r, θ, z), as described above in flowchart 700 in FIG. 7A. Flowchart 700 also describes how the reflectivity M(r, θ, z, a) is found from the recorded signal data us,R(t) for different pairs (S,R) of source transducers and receiver transducers, as a function of time t, measured relative to the starting time of transmission from the source receiver. At 1304, after the correct velocity model has been found, the correct reflectivity M(r, θ, z, a) is found as a function of position (r, θ, z) and reflection angle a. For a given position (r, θ, z), the a dependence of the correct reflectivity is expected to be proportional to:

Here R(a) is the reflection coefficient, a is the reflection angle, C is a background constant that is taken to be 2/3, Rv_P is the reflectivity of the P-wave velocity Vp, Rvs is the reflectivity of the S- wave velocity Vs, and R_P is the reflectivity of the density. For each angle gather (θ, z), the quantities Vs, Vp, and p are defined on each grid-point n of the angle gather. The reflectivity RQ of a quantity Q is defined on each grid-point n as:

At 1306, the α ddpendence of the measured reflectivity M(r, θ, z, a) at a given point (r, θ, z) is fitted to the above expression for R(α), with Rv_P, Rvs, and Rv_s acting as free parameters, using a least squares fit, for example. The best fit values of Rv_P, Rvs, and Rvs provide measured values of the P-wave velocity reflectivity, the S-wave velocity reflectivity, and the density reflectivity. Although this procedure does not immediately give the absolute values of the measured P-wave velocity reflectivity, S-wave velocity reflectivity, and density reflectivity, because M(r, θ, z, a) is in arbitrary units, it does give their relative values, as functions of position (r, θ, z) in the field of view. The P-wave velocity reflectivity, S-wave velocity reflectivity, and density reflectivity are referred to as pre-stack attributes, because they can only be determined from the pre-stack data M(r, θ, z, a), including its dependence on a.

At 1308, the relative values of other pre-stack attributes are optionally found from these pre-stack attributes. For example, the Poisson’s ratio reflectivity R_σ is given by (8/9)(R_{V p} - Rvs), the P-wave impedance reflectivity Ri_p is given by Rv_P + R_P, and the S-wave impedance reflectivity Ris is given by Rvs + R_P. M(r, θ, z, a) may be summed over the different values of a to find the “stack” M(r, θ, z), which may be used as the main output image, and is considered a post-stack attribute. This is optionally found at 1309.

The measured reflectivity M(r, θ, z, a) found at 1304, using the method of flowchart 700, is in arbitrary units. The method of flowchart 700 gives the relative values of measured reflectivity M(r, θ, z, a) at different values of r, θ, z, and a, not the absolute measured reflectivity, and the prestack measured reflectivity attributes found from M(r, θ, z, a) at 1306 and 1308 are also only relative values as a function of r, θ, z. At 1310, one or more of these measured reflectivity attributes are calibrated, to find their magnitudes absolutely. There are several possible methods for doing this calibration, and once it is done for one measured reflectivity, the same calibration can be used for all the measured reflectivities. l)If the ultrasound wave amplitude produced by a given voltage applied to a source transducer is determined experimentally and/or by a simulation, and the ratio of received ultrasound wave amplitude to receiver transducer output voltage is determined experimentally and/or by a simulation, then it may be possible to determine the absolute measured reflectivity M(r, 0, z, a) calculated using the method of flowchart 700, for example using a simulation or a calculation that takes into account the grid spacing of the (r, 0, z) coordinate system, the binning { oci} of reflection angle a, the central ultrasound frequency fo, the corresponding ultrasound wavelength as a function of position (r, 0, z) taking into account the velocity model V(r, 0, z), and the shape of the source function W(T) produced by the source transducers.

2) The post-stack reflectivity M(r, θ, z), calculated from the method of flowchart 700, can be calibrated by looking at its value at the boundaries between bones and soft tissue, where the absolute value of the reflectivity may be determined from known values of P-wave impedance in bone and adjacent soft tissue.

3) If an average value of reflectivity M(r, θ, z), for example the mean absolute value or the root mean squared value, is known for that type of tissue, for example from modeling of some kind, or from some kind of measurements, including post mortem measurements that would be difficult to do on a living subject, then that information can be used to calibrate the measured M(r, θ, z) obtained from a subject, by assuming that the average value for that subject is the same as the known average value.

It should be noted that even if the calibration is initially off by 10% or 20%, when the P- wave impedance is found from the P-wave impedance reflectivity using the inversion procedure described below at 1312, which includes a regularization term that requires the P-wave impedance to be not too different from known average values in different types of body tissue, then the error in calibration is likely to become much smaller.

At 1312, post-stack attributes are calculated from the corresponding calibrated measured reflectivities found at 1310, using an inversion procedure and an optimization algorithm. For example, P-wave velocity is calculated from the P-wave velocity reflectivity, S-wave velocity is calculated from the S-wave velocity reflectivity, density is calculated from the density reflectivity, P-wave impedance is calculated from the P-wave impedance reflectivity, S-wave impedance is calculated from the S-wave impedance reflectivity, and Poisson’s ratio is calculated from the Poisson’s ratio reflectivity. Specifically, this is done, in each case, by finding a set of trial values of the quantity being calculated, for each grid point n along a path of constant (θ, z) for a given angle gather, that minimizes a function involving the trial values of the quantity being calculated, the corresponding calibrated measured reflectivity at each grid point n found at 1310, the shape I l l W(T) of the wavelet transmitted by the source transducers, and a regularization term that keeps the quantity being calculated close to a physically expected value for the body tissue. For example, in order to find the P-wave impedance Ip(ri) along a path of constant (θ, z) for a given angle gather from the calibrated measured P-wave impedance reflectivity A’/_/,(rj) that was determined at 1310, we find a trial set of values of Ip for each grid point n of the angle gather that minimizes the following function:

Here, ti is the vertical travel time along the path of constant (θ, z) from the top of the angle gather down to grid point n. The notation 7p(t) is used in this equation as shorthand for 7p(r(t)), the trial set of values Ip evaluated at the r for which the vertical travel time is t. If r is not one of the grid points n at which //?(r) is defined, then //?(r) is to be evaluated by interpolation among neighboring grid points, using any known interpolation scheme, for example linear or quadratic interpolation. The sum over z includes all the grid points in the angle gather except the last one, because t_I+i is not defined for the last grid point. The wavelet W(T) is the shape of the transmitted signal, with T = 0 being the time of transmission from which the recorded received signals are measured, when they are used in the method of flowchart 700 to find M(r, θ, z, a). For example if the recorded signals are measured from the middle of the transmitted signal, then T = 0 is halfway between the first time at which iifrjdiffcrs from zero, and the last time at which W(T) differs from zero. The wavelet is normalized so that The integral over T in the expression above

is over the range for which the wavelet W(T) is non-zero. The quantity in the first term in square brackets is the P-wave

impedance reflectivity as a function of n, that you would expect to measure if the P-wave impedance has the form of the trial values Ip(ri). It is the convolution of

Ipiti)) for the trial values, with the wavelet function w(t). The first term in square brackets is minimized when the trial values for the P-wave impedance are expected to produce a measured P- wave impedance reflectivity that matches the actual measured P-wave impedance reflectivity A7_/;(ri) obtained using the method of flowchart 700. The second term in square brackets is the difference between the trial P-wave impedance value /p(n) and a physically reasonable value for soft body tissue, 2.7 x 10⁶ kg/m²s. This second term is a regularization term, and when the sum of the squares of the first and second terms is minimized, the regularization term prevents the optimized trial function Ip(t) from being very noisy due to over-determination by noise in the measured P-wave impedance reflectivity Ri_p(t). Optionally, if the angle gather is known to include bone or other tissue that has a very different P-wave impedance, then a reasonable value for that tissue, instead of 2.7 x 10⁶ kg/m²s, is used for grid points n that are expected to be in that tissue. The coefficient C in front of the second term is chosen to provide a good trade-off between minimizing the first term and minimizing the second term, in finding the optimized P-wave impedance The optimum

value of C may make the first and second terms in brackets equal, or at least comparable in magnitude, for the set of trial values that minimizes the function, and this value of C can be

found by trial and error. The same method is used to find other quantities from their corresponding reflectivities, in each case using a regularization term that keeps the quantity close to a reasonable value in body tissue, and choosing an appropriate value of C, the coefficient in front of the regularization term, that is expected to work well for that quantity.

In addition to the post-stack attributes that are found from corresponding reflectivities at 1312 of FIG. 13, other post-stack attributes may be found at 414 in FIG. 4, for example directly from the stacked reflectivity M(r, θ, z), or from other reflectivity attributes. For example, the complex trace of the reflectivity is optionally found at 1314 of FIG. 13, a quantity that is potentially useful because the reflectivity in body tissue tends to oscillate in waves as a function of position. At 1316, one or more other post-stack attributes are found, from other mathematical manipulations of pre-stack or post-stack attributes. The mathematical manipulations may include, for example, taking one or more of a time derivative of any order, a spatial gradient or any other spatial derivative of any order, a temporal filter and a spatial filter of the attribute.

As used herein the term “about” refers to ±10%.

The terms "comprises", "comprising", "includes", "including", “having” and their conjugates mean "including but not limited to".

The term “consisting of’ means “including and limited to”.

The term "consisting essentially of" means that the composition, method or structure may include additional ingredients, steps and/or parts, but only if the additional ingredients, steps and/or parts do not materially alter the basic and novel characteristics of the claimed composition, method or structure.

As used herein, the singular form "a", "an" and "the" include plural references unless the context clearly dictates otherwise. For example, the term "a compound" or "at least one compound" may include a plurality of compounds, including mixtures thereof.

Throughout this application, various embodiments of this invention may be presented in a range format. It should be understood that the description in range format is merely for convenience and brevity and should not be construed as an inflexible limitation on the scope of the invention. Accordingly, the description of a range should be considered to have specifically disclosed all the possible subranges as well as individual numerical values within that range. For example, description of a range such as from 1 to 6 should be considered to have specifically disclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numbers within that range, for example, 1, 2, 3, 4, 5, and 6. This applies regardless of the breadth of the range.

Whenever a numerical range is indicated herein, it is meant to include any cited numeral (fractional or integral) within the indicated range. The phrases “ranging/ranges between” a first indicate number and a second indicate number and “ranging/ranges from” a first indicate number “to” a second indicate number are used herein interchangeably and are meant to include the first and second indicated numbers and all the fractional and integral numerals therebetween.

It is appreciated that certain features of the invention, which are, for clarity, described in the context of separate embodiments, may also be provided in combination in a single embodiment. Conversely, various features of the invention, which are, for brevity, described in the context of a single embodiment, may also be provided separately or in any suitable subcombination or as suitable in any other described embodiment of the invention. Certain features described in the context of various embodiments are not to be considered essential features of those embodiments, unless the embodiment is inoperative without those elements.

Although the invention has been described in conjunction with specific embodiments thereof, it is evident that many alternatives, modifications and variations will be apparent to those skilled in the art. Accordingly, it is intended to embrace all such alternatives, modifications and variations that fall within the spirit and broad scope of the appended claims.

It is the intent of the Applicant(s) that all publications, patents and patent applications referred to in this specification are to be incorporated in their entirety by reference into the specification, as if each individual publication, patent or patent application was specifically and individually noted when referenced that it is to be incorporated herein by reference. In addition, citation or identification of any reference in this application shall not be construed as an admission that such reference is available as prior art to the present invention. To the extent that section headings are used, they should not be construed as necessarily limiting. In addition, any priority document(s) of this application is/are hereby incorporated herein by reference in its/their entirety.

Claims

WHAT IS CLAIMED IS:

1. A method of obtaining a medical ultrasound image of a field of view of a subject’s body with a coordinate system comprising grid points, the method comprising: a) placing ultrasound source transducers and ultrasound receiver transducers on the body of the subject, surrounding the field of view; b) recording ultrasound signals each transmitted by one of the source transducers through the body and received by one of the receiver transducer, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) using a velocity model of the field of view to find travel times to each of a plurality of the grid points from each of the source transducers that transmitted signals, and from each of the plurality of grid points to each of the receiver transducers that received signals; d) constructing gathers as a function of grid point and of discrete bins of a parameter related to ray path length, by, for each grid point in the plurality of grid points, and for each of a plurality of bins of the parameter related to ray path length:

2. A method according to claim 1, wherein the source transducers and receiver transducers are arranged on the body of the subject such that, for every point in a volume comprising most of the field of view, for every plane that contains the point, there are one or more transducers located on each side of the plane.

3. A method according to claim 2, wherein for most of the source transducers, for every point in a volume comprising most of the field of view, for every plane that contains the point, there are one or more receiver transducers located on each side of the plane, for which an ultrasound signal transmitted by that source transducer and received by that receiver transducer is one of the ultrasound signals recorded.

4. A method according to any of the preceding claims, wherein the gathers comprise angle gathers, the parameter related to ray path length is reflection angle at the grid point, and selecting recorded signals comprises finding reflection angle at the grid point for each recorded signal, and selecting signals for which the reflection angle at the grid point is in that bin.

5. A method according to claim 4, where finding reflection angle at the grid point for a recorded signal comprises using the velocity model to find a ray path of ultrasound propagating from the source transducer that transmitted the signal to the grid point and from the grid point to the receiver transducer that received the signal.

6. A method according to claim 4 or claim 5, comprising using the angle gathers to find pre-stack attributes of one or more of Poisson’s ratio reflectivity, P wave velocity reflectivity, P wave impedance reflectivity, S wave velocity reflectivity, S wave impedance reflectivity, and density reflectivity.

7. A method according to any of claims 4-6, wherein constructing the gathers comprises constructing an angle gather for at least one grid point for a reflection angle greater than 60 degrees.

8. A method according to any of the preceding claims, wherein correcting the velocity model comprises: a) for the plurality of grid points, finding a value of a trend of the gather with the parameter, for fixed grid point, that would be expected if there were an error in average velocity from an outer surface of the body inward to the grid point; b) using the trend to find a corrected average velocity as a function of grid point; and c) locally inverting the corrected average velocity as a function of grid point to find a corrected interval velocity as a function of grid point.

9. A method according to any of the preceding claims, also comprising constructing a main reflectivity image by stacking the gathers over values of the parameter.

10. A method according to any of the preceding claims, wherein placing the ultrasound source transducers and the ultrasound receiver transducers on the body of the subject comprises wrapping around a part of the body a flexible stretchable frame, adapted to conform to different body shapes, with the transducers mounted on the frame and held in acoustic contact with the body by the frame.

11. A method according to claim 10, also comprising obtaining information, about relative positions, at times when the ultrasound signals are transmitted or received, of the source transducers and the receiver transducers.

12. A method according to any of the preceding claims, wherein the coordinate system comprises an (r, θ, z) coordinate grid, and constructing gathers comprises constructing the gathers with a discrete set of grid coordinates (θ, z) defining the locations of gathers extending into the field of view away from an outer surface of the subject’s body, and wherein there is a range of values of θ and a range of values of z that include at least ten pairs (θ, z) that belong to the discrete set, and a range of values of r such that all of the gathers defined by said pairs (θ, z) extend over said range of r, with r increasing along each of said gathers in a direction going into the field of view further away from a closest outer surface of the subject’s body, and wherein, over said range of r, the average distance between said gathers at a given r is smaller for larger values of r.

13. A method according to any of the preceding claims, wherein the field of view has a tube-shaped outer surface, corresponding to an outer surface of the subject’s body, and wherein the coordinate system comprises an (r, θ, z) coordinate grid, with θ a periodic azimuthal coordinate around the tube, z an axial coordinate along the tube, and r a coordinate going into the field of view from the outer surface, and wherein constructing the gathers comprises constructing the gathers with a discrete set of grid coordinates (θ, z) defining the locations of the gathers extending into the field of view away from the outer surface, and wherein there is a range of values of z that include at least ten pairs (θ, z) that belong to the discrete set, and a range of values of r such that all of the gathers defined by said pairs (θ, z) extend over said range of r, with r increasing along each of said gathers in a direction going into the field of view away from the outer surface, and wherein, for said range of z and said range of r, the field of view extends over a full range of values of 0.

14. A method according to claim 13, wherein correcting the velocity model using the gathers comprises: a) finding a trend in a dependence of the gather on the parameter for fixed r, θ and z, that would be expected if there were an error in the velocity model, and using the trend to find an average velocity as a function of r for each pair (θ, z) that corrects the error; and b) finding a corrected interval velocity as a function of r, θ, and z, at least over said ranges of r, θ, and z, by locally inverting the average velocity as a function of r for each pair (θ, z), and smoothing the interval velocity, as a function of r and z, at least over 0.

15. A method according to any of the preceding claims, wherein the parameter comprises offset between source transducer and receiver transducer.

16. A method according to any of the preceding claims, also comprising using the prestack attribute to find a corresponding post-stack attribute of at least a portion of the field of view as a function of grid point.

17. A method according to claim 16, wherein the pre-stack attribute comprises a reflectivity of a quantity, the post-stack attribute corresponding to that pre-stack attribute is the quantity, and finding the corresponding post-stack attribute comprises performing a local inversion procedure using an optimization algorithm.

18. A method according to claim 16 or claim 17, wherein using the pre-stack attribute to find a corresponding post-stack attribute comprises finding one or more of Poisson’s ratio, P wave velocity, P wave impedance, S wave velocity, S wave impedance, and density.

19. A method according to any of the preceding claims, comprising finding one or more post-stack attributes by performing a mathematical manipulation of a pre-stack attribute or of a post-stack attribute.

20. A method according to claim 19, wherein performing the mathematical manipulation comprises taking a complex trace.

21. A method of obtaining at least a first medical ultrasound image of a field of view of a first body part of a first subject, the method comprising: a) wrapping a flexible stretchable frame, with ultrasound source transducers and ultrasound receiver transducers mounted on it, around the first body part, with the ultrasound transducers held by the frame in good acoustic contact with an outer surface of the first body part; b) recording a first set of ultrasound signals, each signal transmitted from one of the source transducers through the first body part and received by one of the receiver transducers, while the transducers are held in good acoustic contact with the outer surface of the first body part, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) determining, from the recorded first set of ultrasound signals, information about relative spatial positions of the source transducers that transmitted one of the first set of signals and the receiver transducers that received one of the first set of signals, at times when the first set of ultrasound signals were transmitted or received; and d) reconstructing the first ultrasound image from the recorded first set of ultrasound signals and the information about the relative spatial positions of the transducers at times when the signals were transmitted or received, after determining the information about the relative spatial positions of the transducers.

22. A method according to claim 21, wherein reconstructing the first ultrasound image from the recorded first set of ultrasound signals uses a coordinate system comprising a plurality of grid points, and comprises: a) using a velocity model of the field of view to find travel times to each grid point from each source transducer that transmitted a signal and from each grid point to each receiver transducer that received a signal; and b) constructing gathers of the recorded signals, using the travel times, as a function at least of grid point; and wherein constructing the gathers is done after determining the information about the relative spatial positions of the transducers.

23. A method according to claim 21 or 22, of obtaining also a second medical ultrasound image of a field of view of a second body part, of the first subject or a different subject, the second body part having a shape, size, or both that differs from the first body part, the method also comprising: a) after recording the first set of ultrasound signals, wrapping the frame around the second body part, with the ultrasound transducers held by the frame in good acoustic contact with an outer surface of the second body part, adjusting the frame if necessary to produce the good acoustic contact; b) recording a second set of ultrasound signals each signal transmitted from one of the source transducers through the second body part and received by one of the receiver transducers, while the transducers are held in good acoustic contact with the outer surface of the second body part, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) determining, from the recorded second set of ultrasound signals, information about relative spatial positions of the source transducers that transmitted one of the second set of signals and the receiver transducers that received one of the second set of signals, at times when the second set of ultrasound signals were transmitted or received; and d) calculating the second ultrasound image from the recorded second set of ultrasound signals and the information about the relative spatial positions of the transducers at times when the signals were transmitted or received, after determining the information about the relative spatial positions of the transducers.

24. A method according to any of claims 21-23, wherein the first body part is a trunk of the first subject, and the transducers are mounted on the frame, and the frame is wrapped around the trunk, such that, along at least one path along the frame encircling the trunk, there is no gap between transducers that is wider than 20% of the length of the path.

25. A method according to claim 24, wherein the field of view comprises at least 50% of at least one axial cross-section of the trunk.

26. A method according to claim 25, wherein the field of view is three dimensional, and includes at least 50% of the volume of one or more of: 1) the lungs and the heart; 2) the stomach and the liver; 3) the small intestine and the colon.

27. A method according to claim 25 or claim 26, of obtaining also a second medical ultrasound image of a field of view of a trunk of a second subject, differing in size, shape, or both from the trunk of the first subject, the method comprising: a) after recording the first set of ultrasound signals, wrapping the frame around the trunk of the second subject, with the ultrasound transducers held by the frame in good acoustic contact with an outer surface of the trunk of the second subject, adjusting the frame if necessary to produce the good acoustic contact; b) recording a second set of ultrasound signals, each signal transmitted from one of the source transducers through the trunk of the second subject and received by one of the receiver transducers, while the transducers are held in good acoustic contact with the outer surface of the trunk of the second subject, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; c) determining, from the recorded second set of ultrasound signals, information about relative spatial positions of the source transducers that transmitted one of the second set of signals and the plurality of receiver transducers that received one of the second set of transducers, at times when the second set of ultrasound signals were transmitted or received; and d) calculating the second ultrasound image from the recorded second set of ultrasound signals and the information about the relative spatial positions of the transducers at times when the signals were transmitted or received, after determining the information about relative spatial positions of the transducers.

28. A method according to any of claims 21-27, wherein determining information about relative spatial positions of the transducers comprises determining spatial positions of the transducers with a standard deviation error of less than 0.5 mm.

29. A method according to any claims 21-28, also comprising providing initial estimates for the relative spatial positions of the source transducers that transmitted one of the signals in the first set of signals and the receiver transducers that received one of the signals in the first set of signals, at the times when the first set of ultrasound signals were transmitted or received, wherein determining information about the relative spatial positions of the transducers comprises using the initial estimates.

30. A method according to claim 29, wherein providing initial estimates for the relative spatial positions of the transducers comprises: a) providing information on locations where the transducers are mounted on the frame; and b) using the information on the locations where the transducers are mounted on the frame, and a model of the first subject’s body, to estimate relative spatial positions of the transducers when the frame is wrapped around the first body part.

31. A method according to claim 29 or claim 30, wherein the frame comprises one or more sensors that generate data of a distance between the transducers in one or more pairs of transducers mounted adjacent to each other on the frame, one or more sensors that generate data of a curvature or angle between the transducers in one or more sets of three transducers mounted adjacent to each other on the frame, or both.

32. A method according to claim 31, wherein providing the initial estimates for the relative spatial positions of the transducers comprises using the data generated by said sensors.

33. A method according to claim 32, wherein providing the initial estimates for the relative spatial positions of the transducers also comprises: a) providing information on locations where the transducers are mounted on the frame; and b) using the information on the locations where the transducers are mounted on the frame, and a model of the first subject’s body, to estimate relative spatial positions of the transducers when the frame is wrapped around the first body part.

34. A method according to any of claims 21-33, wherein determining information about the relative spatial positions of the transducers from the recorded first set of ultrasound signals comprises determining information about the relative spatial positions of the transducers from only a sub-set of the recorded first set of ultrasound signals where the source transducer and the receiver transducer are located relatively close to each other, compared to distances between the source transducer and the receiver transducer for ultrasound signals in the recorded first set that are not in the sub-set.

35. A method according to any of claims 21-34, wherein determining information about relative spatial positions of the transducers comprises using a computer running a program configured to improve its performance by machine learning.

36. A method according to any of claims 21-35, wherein determining information about relative spatial positions of the transducers comprises using a computer running a neural network configured to determine the positions from the recorded first set of ultrasound signals.

37. A method according to claim 36, wherein the neural network comprises a convolutional neural network (CNN).

38. A method according to claim 37, wherein the CNN has a U-shaped architecture.

39. A method according to claim 38, wherein running the CNN comprises: a) receiving a primary input data volume comprising the recorded signals from the first set of source transducer-receiver transducer pairs; b) passing the primary input data volume through the U-shaped architecture to generate a first intermediate data volume; c) receiving a secondary input data volume comprising initial estimates of the positions of the source transducers and the receiver transducers in the first set of pairs; d) performing a first combining of data volumes using the first intermediate data volume and the secondary input data volume, by one of:

40. A method according to claim 39, wherein using the second intermediate data volume to generate the output data volume comprises performing a second combining of data volumes using the second intermediate data volume and one of the first intermediate data volume and the secondary input data volume, by one of:

41. A method according to claim 40, wherein using the third intermediate data volume to generate the output data volume comprises performing a third combining of data volumes using the third intermediate data volume and a different one of the first intermediate data volume and the secondary input data volumes than used in the second combining of data volumes, by one of:

42. A method according to claim 41, wherein the first combining of data volumes, the second combining of data volumes, and the third combining of data volumes each use a different one of element-by-element addition, element-by-element multiplication, and concatenation.

43. A method according to any of claims 39-42, wherein each of the source transducers transmits signals to the receiver transducers at a different time, and the output data volume comprises a position of each of the receiver transducers at each of the different times when one of the source transducers transmits signals to it.

44. A method of training a neural network for use in claim 36, or for training a program for use in claim 35, to determine the relative spatial positions of the transducers from the recorded first set of ultrasound signals, the method comprising using numerical solutions of a 3-D wave equation for ultrasound waves transmitted through a body part from a source transducer to a receiver ultrasound transducer to simulate the recorded first set of ultrasound signals, using at least one model of ultrasound propagation properties as a function of position in the body part, for a plurality of different combinations of positions of source transducers and positions of receiver transducers on the outer surface of the body part.

45. A method according to claim 44, wherein using at least one model of ultrasound propagation properties as a function of position in the body part comprises using a plurality of models of ultrasound propagation properties as a function of position in the body part and shape and size of the body part, corresponding to a range of subjects in a population of subjects.

46. A method of obtaining one or more types of medical ultrasound images of bodies of subjects belonging to a population, the method comprising: a) providing a subject belonging to the population, general information about the subject, and a type of ultrasound image to be obtained of a field of view of the subject’s body; b) placing on the subject’s body an array comprising ultrasound source transducers and ultrasound receiver transducers, the transducers in good acoustic contact with an outer surface of the subject’s body, and positioned so that the receiver transducers can receive ultrasound signals transmitted into the field of view by the source transducers; c) automatically providing an initial velocity model for sound speed as a function of position in at least a part of the subject’s body including the field of view, one or more parameters for transmitting a set of ultrasound signals through the subject’s body from a plurality of the source transducers to a plurality of the receiver transducers, one or more parameters for finding travel times using a velocity model, one or more parameters for constructing gathers from recorded received ultrasound signals, one or more parameters for correcting the velocity model, one or more parameters for iterating the correction of the velocity model, one or more parameters for constructing a pre-stack attribute from the gathers, and one or more parameters for reconstructing an ultrasound image from the pre-stack attribute; d) recording a set of ultrasound signals, each signal transmitted by one of the source transducers through the subject’s body and received by one of the receiver transducers, according to the parameters for transmitting the set of ultrasound signals through the subject’s body, associating each recorded signal with an identifier of its source transducer and an identifier of its receiver transducer; e) using a velocity model of the field of view to find travel times to each of a plurality of grid points from each source transducer that transmitted one of the signals and from each of the plurality of grid points to each receiver transducer that received one of the signals, according to the parameters for using a velocity model to find travel times; f) constructing gathers as a function of grid point and of discrete bins of a parameter related to ray path length, by, for each grid point in the plurality of grid points, and for each of a plurality of bins of the parameter related to ray path length:

3) setting the gathers at that grid point and at that bin of the parameter, to a weighted or unweighted sum of the evaluated selected signals; according to the parameters for constructing gathers; g) correcting the velocity model using information, from dependence of the gathers at a given grid point on value of the parameter, about how ray path length affects travel time, according to the parameters for correcting the velocity model; h) correcting the travel times using the corrected velocity model, according to the parameters for using a velocity model to find travel times; i) repeating (f), using the corrected travel times; j) iterating (g), (h) and (i) until a convergence criterion is satisfied, according to the parameters for iterating the correction of the velocity model and the gathers; k) constructing a pre-stack attribute of the field of view as a function of grid point, using information from dependence of the gathers on grid point and on the parameter, according to the parameters for constructing a pre-stack attribute; and 1) reconstructing an image of the field of view from the pre- stack attribute as a function of grid point, according to the parameters for reconstructing an ultrasound image from the pre- stack attribute.

47. A method according to claim 46, wherein the array comprises a flexible stretchable frame on which the transducers are mounted, and placing the array on the body comprises placing the flexible stretchable frame around the body.

48. A method according to claim 46 or claim 47, wherein one or more of automatically providing the initial velocity model, automatically providing one or more parameters for transmitting a set of ultrasound signals through the subject’s body from a plurality of the source transducers to a plurality of the receiver transducers, one or more parameters for finding travel times using a velocity model, one or more parameters for constructing gathers from recorded received ultrasound signals, one or more parameters for correcting the velocity model, one or more paramaters for iterating the correction of the velocity model, one or more parameters for constructing a pre-stack attribute from the gathers, and one or more parameters for reconstructing an ultrasound image from the pre-stack attribute, comprises automatically providing based on one or more of general information about the subject, the type of ultrasound image, and data about spacing of the transducers obtained automatically after placing the frame around the subject’s body.

49. A method according to any of claims 46-48, also comprising automatically providing one or more parameters for finding reflection angle, wherein the parameter related to ray path length comprises reflection angle, and selecting recorded signals comprises finding reflection angle at the grid point for each recorded signal, and selecting signals for which the reflection angle at the grid point is in that bin.

50. A method according to claim 49, also comprising automatically providing one or more parameters for using a velocity model to find ray paths of propagating ultrasound, where finding reflection angle, at the grid point, of an ultrasound signal travelling from a source transducer to the grid point to a receiver transducer comprises using the velocity model to find a ray path of ultrasound propagating from the source transducer to the grid point and from the grid point to the receiver transducer, according to the parameters for using the velocity model to find ray paths.

51. A method according to any of claims 46-50, wherein constructing gathers comprises constructing gathers on an (r, θ, z) coordinate grid, with a discrete set of grid coordinates (θ, z) defining the locations of gathers extending into the field of view away from an outer surface of the subject’s body, and wherein there is a range of values of θ and a range of values of z that include at least ten pairs (θ, z) that belong to the discrete set, and a range of values of r such that all of the angle gathers defined by said pairs (θ, z) extend over said range of r, with r increasing along each of said angle gathers in a direction going into the field of view away from a closest outer surface of the subject’s body, and wherein, over said range of r, the average distance between said angle gathers at a given r is smaller for larger values of r.

52. A method according to claim 51, wherein correcting the velocity model using the gathers comprises: a) finding a trend in a dependence of the gather on the parameter related to path length, for fixed r, θ and z, that would be expected if there were an error in the velocity model, and using the trend to find an average velocity as a function of r for each pair (θ, z) that corrects the error; and b) finding a corrected interval velocity as a function of r, θ, and z, at least over said ranges of r, θ, and z, by locally inverting the average velocity as a function of r for each pair (θ, z), and smoothing the interval velocity, as a function of r and z, at least over 0.

53. A system for obtaining a medical ultrasound image of a field of view of a body part of a subject with a coordinate grid comprising grid points, the system comprising: a) an ultrasound array comprising:

1) a flexible stretchable frame of a two-dimensional structure;

3) a plurality of ultrasound receiver transducers, one or more of which are also source ultrasound transducers or none of which are also source ultrasound transducers, mounted on the frame in a distributed arrangement over an area of the two-dimensional structure; and

4) a fastener to hold different ends of the two-dimensional structure of the frame together, holding the frame in place when the frame is wrapped around the body part; wherein the transducers are mounted on the frame with faces facing a same direction normal to the area of the two-dimensional structure when the frame is unwrapped and flat, such that the transducers are held pressed in good acoustic contact with an outer surface of the body part when the frame is wrapped and fastened around the body part in at least one configuration, the array adapted for wrapping and fastening around the body part of the subject with the transducers held pressed in good acoustic contact with an outer surface of the body part; b) a controller, comprising an output interface, configured to:

6) correct the travel times using the corrected velocity model;

7) repeat (4), using the corrected travel times;

8) iterate (5), (6) and (7) until a convergence criterion is satisfied;

9) construct a pre-stack attribute of the field of view as a function of grid point, using information from dependence of the gathers on grid point and on the parameter;

54. A system according to claim 53, wherein the ultrasound array also comprises sensors incorporated in the frame that measure distance along the frame between one or more pairs of the transducers, angles in the two-dimensional structure of the frame between one or more sets of three of the transducers, or both, when the frame is wrapped around and the body part and fastened with the transducers held pressed against the outer surface of the body part.

55. A method of obtaining a medical ultrasound image of a field of view of a subject’s body with grid points, comprising: a) placing ultrasound source and receiver transducers on the body, surrounding the field of view; b) recording ultrasound signals transmitted by a source transducer through the body and received by a receiver transducer; c) using a velocity model to find travel times between grid points and source and receiver transducers; d) constructing gathers as a function of grid point and a parameter related to ray path length, using pre-stack depth Kirchhoff migration; e) correcting the velocity model using information from the gathers about how ray path length affects travel time; f) constructing a pre- stack attribute of the field of view as a function of grid point, from the gathers; and g) reconstructing an image of the field of view from the pre- stack attribute as a function of grid point.