The present invention relates to an improved method of, and apparatus for, non-invasive imaging of regions of the body using ultrasound-generated data. More particularly, the present invention relates to a method for non-invasive imaging of bone- or gas-containing regions of the body using waveform inversion of ultrasound generated data, which may be useful in the diagnosis and prognosis of pathologies, particularly brain pathologies such as stroke.
Medical ultrasound is a well established technology that uses high frequency acoustic waves (generally at frequencies in excess of 20 kHz) to generate acoustic images of the human or animal body. The application of this technique ranges from obstetric ultrasound to identification of pathologies, echocardiography or to drive interventions in real time.
Different types of images can be generated using ultrasound technology. 2D acoustic impedance maps can be generated. Alternatively, blood flow maps or motion of tissue over time can be measured and recorded using Doppler effects.
The advantages of ultrasound are many—in particular, this technique provides a real time, portable and cheap imaging technology which does not require ionising radiation. In the majority of medical applications, reflected acoustic energy is detected to generate acoustic impedance maps. Reflection mode imaging such as this is very successful in particular areas of medical imaging—for example, the imaging of a foetus in the womb, or the detection of breast cancer. The relative simplicity of such imaging enables images to be provided in real-time because the algorithms that generate such images are only based on a delay and stack strategy. The only information used to create the image is the time at which the echoes arrive at the receiver, and the strength of those echoes.
However, such techniques have met with limited success when attempting to image regions of the body comprising bone or gas pockets. This is because such materials within the body cause scattering, attenuation and phase distortion of ultrasound signals, making imaging difficult. In addition, imaging of the brain is particularly difficult due to the presence of a strong ultrasound reflector (the skull) in close proximity to the ultrasound source.
For example, “Ultrasound reflection mode computed tomography through a skullbone”, J Ylitalo, IEEE transactions on biomedical engineering Vol 37, No. 11 (November 1990) discloses an ultrasonic reflection mode CT method which showed some success in respect of paediatric brain diagnosis however was unreliable when applied to thicker adult skulls.
An alternative technique to reflection mode analysis is transmission mode analysis. “Computerized ultrasound tomography of the human head: Experimental results”, K. A. Dines et al, Ultrasonic Imaging, 3, 342-351 (1981) discloses methods for using transmission analysis human skulls. However, again, whilst encouraging results were seen for paediatric skulls, the technique proved unreliable when applied to adult skulls.
A further alternative ultrasound technique is to utilise dispersion to categorise brain injuries. U.S. Pat. No. 8,834,376 discloses methods for utilising dispersion to provide information on the composition of inter-cranial tissues. However, such a method cannot provide imaging of brain regions.
To date, the only successful techniques for imaging the brain using ultrasound have relied upon obtaining imaging through thinner regions of the skull, either using naturally-thinner skulls (e.g. a young child's skull) or through natural or man-made (e.g. bored) acoustic windows within the skull. The same applies to other areas of the body containing significant sources of acoustic absorption, reflection or dispersion, or heterogenous regions (e.g. bone or gas regions or interfaces therebetween). Therefore, there exists a technical problem in the art that ultrasound imaging of such regions cannot reliably be performed using known techniques.
According to a first aspect of the present invention there is provided 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, the method comprising the steps of: a) providing an ultrasound observed data set derived from a measurement of one or more ultrasound waveforms generated by at least one source of ultrasound energy, the ultrasound energy being detected by a plurality of receivers located at an opposing side of a region within the intra-cranial cavity with respect to at least one source such that the receivers detect ultrasound waveforms from the source which have been transmitted through the skull and intra-cranial cavity, the observed data set comprising a plurality of observed data values; b) providing at least one starting model for at least a portion of the head comprising a skull component and a soft tissue component, the skull component comprising a plurality of model parameters representative of the physical properties and morphology of the skull through which intra-cranial tissue is being imaged, and the soft tissue component comprising a plurality of parameters representative of the physical properties of the intra-cranial tissue being imaged; c) generating a predicted data set comprising a plurality of predicted data values from the starting model of the skull and of the intra-cranial tissue; d) comparing the observed and predicted data values in order to generate an updated model of at least one physical property within at least a region of the intra-cranial cavity; and e) using the updated model to image a region of the inter-cranial cavity to identify tissue composition and/or morphology within the intra-cranial cavity.
In one embodiment, step b) comprises: f) acquiring subject data relating to the subject, and providing at least the skull component of the starting model based on the acquired subject data.
In one embodiment, the acquired subject data is obtained from a measurement performed on the subject and/or from empirical data relating to the subject.
In one embodiment, the skull component is selected from a group of predetermined skull components based on the acquired subject data.
In one embodiment, the skull component is selected from a group of predetermined skull components based at least in part upon a matching process between at least a part of the observed data set and a group of starting predicted data sets generated from the respective group of the skull components of the starting models.
In one embodiment, one or more skull components of the starting model are generated from measured experimental data.
In one embodiment, the skull component of the starting model is generated based on experimental data from one or more of the following: reflection ultrasound; low-frequency transmitted ultrasound; X-ray computed tomography; shear sensors attached to the head of a subject; laser measurement of the head of a subject; and physical measurement of the head of the subject.
In one embodiment, step b) further comprises: g) processing at least a part of the observed data set to generate and/or refine at least the skull component of the starting model.
In one embodiment, the ultrasound data set is derived from a measurement of ultrasound waveforms generated by a plurality of sources of ultrasound energy, the ultrasound energy being detected by a plurality of receivers, wherein the sources and receivers are located such that the receivers detect transmitted ultrasound waveforms from the sources which have been transmitted through the skull and inter-cranial cavity and/or reflected ultrasound waveforms that have been reflected by the inner and/or outer boundaries of the skull.
In one embodiment, at least the reflected waveforms of the observed data set are used to recover a numerical model of the geometry of at least a part of the skull, at least a part of the skull component of the starting model provided in step b) being derived from the numerical model.
In one embodiment, analysing at least the transmitted waveforms of the said observed dataset in order to recover a numerical model of at least one physical property within at least a region of the intra-cranial cavity, and analysing both reflected and transmitted waveforms in order to recover at least one physical property of the skull, by comparison of the observed reflected and transmitted waveforms with predicted waveforms that have been simulated numerically and/or generated experimental using at least one numerical and/or physical and/or in vivo predicted model for which the relevant geometry and property or properties are known and/or can be inferred or approximated.
In one embodiment, the observed data set comprises a plurality of measurements of one or more ultrasound waveforms generated by at least one source of ultrasound energy, wherein each measurement is taken in a plane.
In one embodiment, one or more planes intersect.
In one embodiment, one or more planes are parallel and offset with respect to each other.
In one embodiment, at least a portion of the said observed and predicted waveforms differ in phase by more than half a cycle at the lowest frequency present in the said observed dataset.
In one embodiment, one or more ultrasound sources emit ultrasound energy having one or more frequencies in the region of 50 kHz to 5 MHz.
In one embodiment, one or more ultrasound sources emit ultrasound energy having a finite bandwidth.
In one embodiment, step d) is performed using full waveform inversion analysis.
In one embodiment, the skull component of the starting model comprises elements having an acoustic velocity in excess of 2300 m/s.
In one embodiment, the soft tissue component of the starting model comprises elements having an acoustic velocity within the range of 700 to 2300 m/s.
In one embodiment, the soft tissue component of the starting model comprises elements having an acoustic velocity within the range of 1400-1750 m/s.
According to a second aspect of the present invention, there is provided a non-invasive method of generating image data of a body part of a subject using ultrasound energy that is transmitted through the body part of the subject, the body part containing at least one interface between bone, soft tissue and/or gas, the method comprising the steps of: a) providing an ultrasound observed data set derived from a measurement of one or more ultrasound waveforms generated by at least one source of ultrasound energy, the ultrasound energy being detected by a plurality of receivers located at an opposing side of a region within the body part with respect to at least one source such that the receivers detect ultrasound waveforms from the source which have been transmitted through the body part, the observed data set comprising a plurality of observed data values; b) providing at least one starting model representative of the body part being imaged, the starting model comprising first and second components, the first component comprising a plurality of model parameters representative of the physical properties and morphology of the bone and/or gas within the body part of the subject to be imaged and having at least one modelled region having an acoustic velocity below 700 m/s and/or above 2300 m/s and the second component comprising a plurality of parameters representative of the physical properties of the soft tissue within the body part of the subject to be imaged; c) generating a predicted data set comprising a plurality of predicted data values from the starting model; d) comparing the observed and predicted data values in order to generate an updated model of at least one physical property within at least a region of the body part; and e) using the updated model to image a region of the body part to identify tissue composition and/or morphology within the body part.
In one embodiment, step b) comprises: f) acquiring subject data relating to the subject, and providing at least the first component of the starting model based on the acquired subject data.
In one embodiment, the acquired subject data is obtained from a measurement performed on the subject and/or from empirical data relating to the subject.
In one embodiment, the first component is selected from a group of predetermined components based on the acquired subject data.
In one embodiment, the first component is selected from a group of predetermined first components based at least in part upon a matching process between at least a part of the observed data set and a group of starting predicted data sets generated from the respective group of the first components of the starting models.
In one embodiment, one or more first components of the starting model are generated from measured experimental data.
In one embodiment, the first component of the starting model is generated based on experimental data from one or more of the following: reflection ultrasound; low-frequency transmitted ultrasound; X-ray computed tomography; shear sensors attached to the body part of a subject; and physical measurement of the body part of the subject.
In one embodiment, step b) further comprises: g) processing at least a part of the observed data set to generate and/or refine at least the first component of the starting model.
In one embodiment, the ultrasound data set is derived from a measurement of ultrasound waveforms generated by a plurality of sources of ultrasound energy, the ultrasound energy being detected by a plurality of receivers, wherein the sources and receivers are located such that the receivers detect transmitted ultrasound waveforms from the sources which have been transmitted through the body part and/or reflected ultrasound waveforms that have been reflected by any inner and/or outer boundaries of the body part.
In one embodiment, at least the reflected waveforms of the observed data set are used to recover a numerical model of the geometry of at least a part of the body part, at least a part of the first component of the starting model provided in step b) being derived from the numerical model.
In one embodiment, analysing at least the transmitted waveforms of the said observed dataset in order to recover a numerical model of at least one physical property within at least a region of the body part, and analysing both reflected and transmitted waveforms in order to recover at least one physical property of the body part, by comparison of the observed reflected and transmitted waveforms with predicted waveforms that have been simulated numerically and/or generated experimental using at least one numerical and/or physical and/or in vivo predicted model for which the relevant geometry and property or properties are known and/or can be inferred or approximated.
In one embodiment, the observed data set comprises a plurality of measurements of one or more ultrasound waveforms generated by at least one source of ultrasound energy, wherein each measurement is taken in a plane.
In one embodiment, one or more planes intersect.
In one embodiment, one or more planes are parallel and offset with respect to each other.
In one embodiment, at least a portion of the said observed and predicted waveforms differ in phase by more than half a cycle at the lowest frequency present in the said observed dataset.
In one embodiment, one or more ultrasound sources emit ultrasound energy having one or more frequencies in the region of 50 kHz to 5 MHz.
In one embodiment, the one or more ultrasound sources emit ultrasound energy having a finite bandwidth.
In one embodiment, step d) is performed using full waveform inversion analysis.
In one embodiment, the starting model in step b) is at least partly derived from X-ray CT measurement.
In one embodiment, the X-ray CT measurement is performed on the body part to be imaged.
In one embodiment, the X-ray CT measurement and ultrasound measurement are performed simultaneously or sequentially on the subject.
In one embodiment, step a) further comprises: h) utilising at least one source of ultrasound energy to generate one or more ultrasound waveforms; and i) performing a measurement of said one or more ultrasound waveforms utilising a plurality of receivers located at an opposing side of a region within the intra-cranial cavity with respect to the at least one source such that the receivers detect ultrasound waveforms from the source which have been transmitted through the skull and intra-cranial cavity.
In one embodiment, step a) further comprises: j) generating an observed data set from the measurement in step i).
In one embodiment, step a) further comprises: h) utilising at least one source of ultrasound energy to generate one or more ultrasound waveforms; and i) performing a measurement of said one or more ultrasound waveforms utilising a plurality of receivers located at an opposing side of a region within the body part with respect to the at least one source such that the receivers detect ultrasound waveforms from the source which have been transmitted through the body part.
In one embodiment, step a) further comprises: j) generating an observed data set from the measurement in step i).
According to a third aspect of the present invention, there is provided a computer system comprising a processing device configured to perform the method of the first or second aspects.
According to a fourth aspect of the present invention, there is provided a computer readable medium comprising instructions configured when executed to perform the method of the first or second aspects.
According to a fifth aspect of the present invention, there is provided a computer system comprising: a processing device, a storage device and a computer readable medium of the third aspect.
Embodiments of the present invention will now be described in detail with reference to the accompanying drawings, in which:
The present invention, in embodiments, relates to a novel method for imaging of body structures using ultrasound. In embodiments of the present invention, both transmitted and reflected energy is recorded. Using both these techniques, it is possible to obtain information relating to tissues and body structures that lie behind bone (and/or gas) and generate images of such.
The method of the present invention uses Full Waveform Inversion (FWI) or a variant thereof.
Whilst FWI has limitations with regard to providing real-time images, developments in the implementation and processing of the images should bring down the turn-around time by at least an order of magnitude in due course. This is offset by the potential advantages of the high resolution of the images generated by FWI, which may potentially rival MRI.
The experimental configuration 10 comprises a ring 12 which includes at least one ultrasound source 14 and a plurality of receivers 16. The ring 12 is arranged around a head 18 of a subject. The source 14 and/or receivers 16 are acoustically coupled to the head 18. This may be achieved by either locating the source 14 and/or receivers 16 directly on the surface of the head. Alternatively, the head 18 may be immersed in a suitable acoustic medium (such as water) to enable the required coupling.
As shown in
For example, in one embodiment, the ring 12 may translate in one or more directions to obtain data in a number of parallel, offset planes to capture “slices” through the head.
In addition, whilst a ring 12 is shown in
The source 14 generates acoustic ultrasound waves having sufficient vibrational energy to penetrate the skull of the subject and generate sufficient return signals to aid useful detection. The source 14 may be any suitable ultrasound generator. The skilled person will be readily aware of the type of generators that are suitable for use with the acquisition approach of the present invention.
The source 14 is, in embodiments, a point source or a close approximation of a point source. In other words, the source 14 is a source of ultrasonic waves which emits in all directions equally. In the context of the experimental configuration 10, this may be understood to be an isotropic source for the purpose of the detectors 16 such that the source 14 is an isotropic radiator in directions of interest for the measurement (e.g. within the ring 12). It is not essential or necessary that the source 14 is isotropic, or even generates a signal at all, in directions which do not intersect the target structure (e.g. directions away from the head).
However, whilst point sources (or functionally similar sources) are desirable, this need not be the case. Instead, the source 14 may have some directionality within the target region. Whilst it is generally preferred to use isotropic or quasi-isotropic emitters, there may be situations where more focused beams (for example, in particular planes or at particular distances) may be useful. However, in contrast to known ultrasound imaging arrangements, the distortion of signals is a useful parameter which contains information that can be recovered and used to provide imaging data. In contrast, known imaging arrangements are configured to avoid distortion due to the difficulty in correcting for such effects.
As shown in
Whilst a single source 14 and multiple receivers 16 are shown in
In use, ultrasonic waves generated by the source 14 propagate into the head 18 of the subject. The waves are transmitted and refracted through the layers of bone and/or brain matter and/or reflected off the interfaces between them and/or scattered from other heterogeneities within the head and a plurality of return signals is detected by the detectors 16.
An observed data set comprises a measurement, by the multiplicity of receivers 16, of transmitted, reflected and/or refracted acoustic waves originating from the source 14. In general, a partial reflection of the acoustic wave occurs at a boundary or interface between two dissimilar materials, or when the elastic properties of a material changes. Each detected return signal forming an ultrasound trace has an approximate travel time from the source 14 which, for reflected waves, is a two-way travel time from the source 14 to the reflecting element (for example, the interior surface of the skull opposite from the source 14) and back to the respective detector 16.
In the experimental configuration 10 shown in
The time-variation of the reflected and transmitted waves (i.e. the waveforms detected by the receivers 16) is recorded during a pre-determined time period. This time period is selected to be sufficient to capture both reflected and transmitted arrivals that contain information of the properties of the head. A typical value may be of the order of 300-500 μs.
The source 14 may emit ultrasonic waves at any desired frequency, plurality of discrete frequencies or continuous band of frequencies (e.g. a broadband signal). In general, a multiplicity of frequencies is used to provide greater resolution and detail in the produced images. Whilst high frequencies provide greater resolution of smaller physical features, the penetration depth of lower frequencies is better for imaging further through the skull. The present invention is operable to use a range of frequencies which can be resolved together or individually as required.
In embodiments, the frequencies span a range from 400 kHz to 1.3 MHz. However, other ranges could be used with the present invention, and a continuous bandwidth and/or discrete frequency selections ranging from approximately 50 kHz up to 5 MHz could be used with the present invention.
Some or all of the multiple sources 14 (if present) may be activated simultaneously to generate a single large source gather. Alternatively, the sources 14 may be actuated individually.
An observed ultrasonic data set is then acquired by recording the waveforms at the receivers 16 after emission by the one or more sources 14. The observed data set may comprise a plurality of waveform traces. For example, there will be a single waveform trace for each source/receiver combination. These traces form the observed ultrasonic data set.
A predicted data set is obtained by modelling data from a starting model of the acoustic properties of the area investigated, such as wave speed velocity. 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 minimising the least-squares norm between observed and predicted data.
Key to this process of modelling and imaging a region of a subject is the ultrasonic velocity Vp. In a portion of the volume of a subject, V p may be estimated in various ways.
The observed ultrasonic data set is then used as part of a waveform inversion process to extract acoustic properties of the tissues forming the patient's head. An example of full waveform inversion (FWI) will now be described.
FWI is a known method for analysing data, particularly in the field of seismology. FWI is able to produce models of physical properties such as V p in the measured region that have high fidelity and that are well resolved spatially. FWI seeks to extract the acoustic properties of the imaged region of the head from the recorded observed data set. A detailed velocity estimate can be produced using an accurate model with variations on the scale of the ultrasonic wavelength.
The FWI technique involves generating a two or three dimensional model to represent the measured portion of the subject's head or body region and attempting to modify the properties and parameters of the model to generate predicted data that matches the experimentally obtained ultrasonic trace data.
The predicted data is calculated from the model typically using the full two-way wave equation. FWI is an iterative process requiring a starting model. A sufficiently accurate starting model for FWI may be provided by travel-time tomography.
FWI can extract many physical properties (Vp and shear-wave velocities, attenuation, density, anisotropy) of the modelled portion of the subject's body. However, Vp, the P-wave velocity, is a particularly important parameter which the subsequent construction of the other parameters depends heavily upon. Nevertheless, other parameters may be used with the present invention, either alone or in combination. Attenuation and density may also be important parameters in the context of medical imaging. The nature and number of parameters used in a model of a portion of the subject's body will be readily apparent to the skilled person.
The FWI technique seeks to obtain an accurate and high resolution model of the measured region of the subject's body which generates predicted data that matches the recorded data. As set out above, determination of V p is a focus of the technique. However, other parameters such as density and attenuation may also be modelled. Predicted data is calculated using the full two-way wave equation. This is known as the forward problem. This equation can be in the time domain, the frequency domain, or other suitable domains, and it may be elastic or acoustic, isotropic or anisotropic, and may include other physical effects such as attenuation and dispersion. In most cases FWI proceeds using the acoustic approximation with a single component modelled wavefield.
An example of the FWI process in accordance with an embodiment of the present invention will now be described. To test the applicability of the method of the present invention, a target model was developed. The target model is representative of a model of a human head and is used as a target for the FWI process. In other words, the target model is used to generate an example synthetic observed data set, and from this, a starting model is iteratively modified using the FWI process to arrive at a final model. The final model is then compared to the target model to validate the process.
The target model is a 2D synthetic model of a human head and is shown in
Note that the skull has velocities beyond the range of the scale of
A single source gather from a synthetic observed data set generated from the target model of
In general, the parameters of the model are estimated at a plurality of points set out in a grid or volume, but they may be estimated from any suitable parameterisation. The model is used to generate a modelled representation of the ultrasound data set, known as the predicted data set. The predicted data set is then compared to the real-world experimentally obtained observed data set. Then, through use of numerical iteration, the parameters of the model are modified until the predicted data set generated from the model matches the actual observed data to a sufficient degree of accuracy or until sufficient convergence is obtained. This will be explained below.
A general method to update a model will now be described. FWI typically operates on the principle of iteratively updating the starting model to minimise or maximise an objective function through repeated steepest-descent direction calculation, or an analogous technique. An objective function represents some measure of the mismatch or some measure of similarity between the recorded data and the predicted data. A measure of mismatch, obtained for example by subtracting two traces, should be minimised; whereas a measure of similarity, obtained for example by cross-correlating two traces, should be maximised.
Due to the non-linearity in the relationship between the model and the data, the objective function used in FWI will oscillate with changes in the model. This makes it necessary to have a sufficiently accurate starting model for global minimum convergence. The objective function can be formulated in the frequency domain, the time domain or other suitable domain. The choice of domain allows the use of pre-conditioning on either the data or the model update direction that could improve convergence or reduce the non-linearity of the inverse problem.
Frequency domain inversion is equivalent to time domain inversion if all the frequencies are inverted simultaneously. However, the global minimum broadens at lower frequencies reducing how accurate the starting model needs to be for local optimisation method to be successful.
A starting model requires at least two components: (1) A component of the model that represents regions, for example bone or gas, that have values of Vp that differ substantially from values that are typical of soft tissue within the body, and (2) a component of the model that represents values typical of soft tissue within the body. The component of the body represented by component (1) will typically have values of Vp that lie below 700 m/s or above 2300 m/s whereas the portion of the body represented by component (2) will typically have values of Vp that lie close to 1500 m/s. With regard to soft tissue, the acoustic velocity varies from around 1450 m/s for fat to about 1730 m/s for skin. Obtaining a satisfactory starting model for component (2) will normally be straightforward whereas obtaining a satisfactory starting model for component (1) will normally be both essential and more difficult.
An example of a basic starting model for the head of a subject is shown in
The starting model makes no assumptions in relation to the brain's velocities and it is a simple homogeneous model of 1500 m/s. The true velocities for the skull are used in the starting model in this example. The skull corresponds to starting model component (1) and the brain to starting model component (2).
A total of 450 transducers equally spaced around the skull are used and they can act as sources or receivers. Since each source activates one of the transducers while the remaining 449 act as receivers, and all transducers are used as sources resulting in a total of 450 independent experiments.
Commonly, localised gradient-based methods are used with FWI. These methods iteratively update an existing model in a direction that derives from the objective function's direction of steepest descent.
There are numerous ways to quantify the difference (also known as the residual) between the data sets. However, amongst the most common is a least-squares formulation where the sum of the squares of the differences between the two data sets is minimised over all sources and receivers and over all times. In other words, a model is sought that minimises the L2-norm of the data residuals.
The L2-norm expresses the misfit between the two data sets as a single number. This parameter is known as the objective function, although it often takes other names such as the misfit function, the cost function or the functional. The objective function f is a real, positive, scalar quantity, and it is a function of the model m.
In practice, a factor of a half is often included in the definition of the objective function f because it makes later results simpler. The objective function ƒ(m) is shown in equation 1):
where ns, nr and nt are the number of sources, receivers and time samples in the data set.
In the time-domain the data and the data residuals will be real quantities. However, in the frequency domain the data will in general be complex, as will the source and sometimes the model properties. Equation 1) is correctly recited for complex data. However, in the frequency domain, nt would be expressed as nf, i.e. a summation over frequency rather than time.
FWI is a local iterative inversion scheme. FWI comprises numerous different methods. The method described herein is one possible implementation of an FWI method suitable for use with the present invention. However, the skilled person would readily be aware of alternative methods that could be used with the present invention.
A starting model m0 that is assumed to be sufficiently close to the true, ideal model is prepared. The process then makes a series of step-wise improvements to this model which successively reduces the objective function towards zero. Therefore, across an iterative step of the calculation, the objective function needs to be considered for a starting model m0 and a new model m=m0+δm.
For a scalar function of a single scalar variable, the Taylor series can be used, truncated to second order. This generates equation 2):
Differentiating this express with respect to m, and setting the result to zero to minimise f with respect to m=m0+δm, equation 2) becomes:
Then, neglecting second order terms, equation 4) can be derived which expresses the updated to the model δm:
Where ∇mƒ is the gradient of the objective function ƒ with respect to the model parameters, and H is the Hessian matrix of second differentials, both evaluated at m0.
If the model has M parameters, then the gradient is a column vector of length M and the Hessian is an M×M symmetric matrix.
If the number of model parameters M is large, then calculating the Hessian is computationally demanding, and inverting the Hessian exactly is normally computationally intractable. Consequently the method that is typically used is to replace the inverse of the Hessian in equation (9) by a simple scalar a (referred to as the step length). Equation 4) can then be expressed as:
Based on equation 5), conventional FWI can use the method of steepest descent. This involves essentially 5 steps:
5. Iterate from step 2 using the new model until the objective function is minimised.
To calculate the gradient and determine the step length, a Jacobian matrix is used as set out in equation 6):
where J is the Jacobian matrix.
A wave equation for a predicted data set d generated by a source s can be written as:
Ap=s 7)
Where the data set d is a subset of the full wavefield p extracted using the diagonal matrix D that has non-zero unit values only where there are observed data. That is, as set out in equation 8);
d=Dp 8)
Equation 7) can then be differentiated with respect to m, which is equal to zero because s and m are independent:
Equation 9) is then pre-multiplied by the matrix D to extract the wavefield only at the points where data exists. The Jacobian can then be rewritten as:
From equation 10), an expression for the gradient can be derived by recognising that DTδd=δd and substituting equation 10) into equation 6), to derive an expression for the gradient:
So to find the gradient, the forward wavefield p is calculated, the numerical operator A is differentiated with respect to the model parameters and the final term of equation 11) is calculated, which represents a back-propagated residual wavefield.
These terms are then multiplied together for all times and all sources, and summed together to give a value corresponding to each parameter within the model, typically to give one value of the gradient at each grid point within the model.
The final term in equation 11) can be written to arrive at equation 12):
A
T
δp=δd 12)
Equation 12) simply describes a wavefield p that is generated by a (virtual) source Od, and that is propagated by the operator AT which is the adjoint of the operator in the original wave equation. So the term that we need to compute in equation 11) is just the solution of a modified wave equation with the data residuals used as a source.
It is then necessary to compute the step length α. Starting from a current model m0 that generates data d0 and residuals δd0, a new model m1=m0+δm1 that generates data d1 and residuals δd1, where δm is a small change in the opposite direction to the gradient.
Therefore, the aim is to find a new model mα=m0+αδm that generates residuals δmα, where α minimises:
Assuming a linear relationship:
δdα=δd0+α(d1−d0)=δd0+α(δd1−δd0) 14)
By rearranging and differentiating with respect to a, setting the differential equal to zero and solving for α, equations 15) and 16) can be derived:
So, to calculate the step length, a forward calculation is made with a perturbed model, and the residuals from both the original and the perturbed models combined to form equation 15).
Once α has been found, the original starting model m0 can be replaced by mα and the step of the iterative calculation is complete. This process can then be repeated.
Note that iteration is necessary because the problem to be solved is non-linear and the inverse problem has been linearised in particular stages. Implicitly, the method invokes the Born approximation. The Born approximation assumes that the perturbation to a wavefield produced by changing a model is linearly related to the change in the model. This is equivalent to considering only first-order scattering by the perturbation.
Ideally, the above method will lead to a convergence to a model which is a correct representation of the skull of the subject under investigation. However, there are some difficulties associated with obtaining correct convergence.
As set out above, FWI methodology for the objective function above relies upon a gradient decent method to solve the inverse problem. This requires that the starting model should match the observed travel times to within half a cycle. However, real ultrasonic data are limited in their frequency bandwidth. This means that real ultrasonic signals are oscillatory.
An inaccurate starting model may predict data that are more than half a cycle in error with respect to the observed data. Such a situation is described as “cycle skipped”. When this occurs, because the methodology seeks only a local minimum, FWI will tend to modify the model such that the predicted and observed data are brought into alignment at the nearest cycle, and this will neither correspond to their correct alignment nor to the correct model. This misalignment to the nearest cycle will reduce the data misfit, and typical FWI schemes will become stranded in this position—they will have become stuck in a local minimum in the data misfit function rather than being able to find the global minimum which corresponds to the true model.
Other variations of the functional can be combined with the method of the present invention to produce high quality recovered models even when the variations in sound speed exceed in fifty percent of that of the typical soft tissue. Given that different tissues present different acoustic properties, the recovered models of acoustic properties can then be used as diagnostic tools.
Consequently, the method of the present invention has the ability to resolve soft tissue in situations where presence of bone and/or gas would prevent other ultrasound methods from obtaining images of sufficient resolution or detail to be useful in a diagnostic context.
Whilst the above example has been illustrated in relation to an ultrasound scan of the head (i.e. skull and brain), this is not to be taken as limiting. The present invention has applicability to regions of the body which traditionally cannot be imaged using conventional ultrasound methods. In other words, the present invention can be used to image areas of the body containing bone and/or gas which would cause significant difficulties for conventional ultrasound imaging. The present invention is of application in situations where body tissues to be imaged are formed from tissue which has a speed of sound in said material which lies outside the range of 700 to 2300 m/s. The bones of the skull, bones in general, air, gas, metal, and most medical implants lie outside this 50% range of 700 to 2300 m/s with respect to the speed of sound in soft tissue. Note that soft tissue has a speed of sound therein within a few percent of 1500 m/s (around 1450 m/s for fat to about 1730 m/s for skin), whereas air has an approximate acoustic velocity of 350 m/s and bone 3000 m/s.
A generalised method of an embodiment of the invention will now be described with reference to
Initially, it is necessary to obtain a set of experimentally gathered data in order to initiate the imaging procedure. This may be gathered by an experimental arrangement such as the set up shown and described with reference to
Any suitable region of the body of a subject may be imaged using the experimental arrangement. Whilst the example above illustrates the head of a subject being imaged, this need not be the case. The present invention can be used to image areas of the body containing bone and/or gas which would cause significant difficulties for conventional ultrasound imaging, for example, the stomach, arms, legs, liver, chest and other regions not comprised solely of substantially homogeneous soft tissue.
The imaging may be done with any configuration of sensor. If the arrangement of
The gathered ultrasonic data may be optionally pre-processed in various ways including by propagating numerically to regions of the model where experimental data have not been acquired directly. A person skilled in the art would be able to design and undertake such pre-processing as might be necessary or desirable. After such pre-processing, the resultant experimentally gathered ultrasonic data set is known as an “observed ultrasonic data set”.
The observed ultrasonic data set may comprise multiple source 14 emissions. The data comprises pressure as a function of receiver position (on the x-axis) with respect to time (on the y-axis).
The trace data comprises a plurality of observed data points. Each measured discrete data point has a minimum of seven associated location values—three spatial dimensions (x, y and z) for receiver (or detector) position (r), three spatial dimensions (x, y, z) for source location (s), and one temporal dimension measuring the time of observation relative to the time of source initiation, together with pressure magnitude data. The seven coordinates for each discrete data point define its location in space and time. It also comprises one or more measurement parameters which denote the physical property being measured. In this embodiment, a single measurement parameter, pressure is measured. The observed data set is defined as dobs(r,s,t) and, in this embodiment, is in the time domain. For clarity, the following discussion considers a single source-receiver pair and so r, s are not needed.
The actual gathering of the ultrasonic data set is described here for clarity. However, this is not to be taken as limiting and the gathering of the data may or may not form part of the present invention. The present invention simply requires a real-world observed ultrasonic data set upon which analysis can be performed to facilitate medical imaging of a bone- or gas-containing region of the body of a subject.
The method now proceeds to step 102.
Step 102: Provide starting model
At step 102, an initial starting model of the region of the body to be imaged is provided. The model may be provided in either a two dimensional or a three dimensional form. Whilst the illustrated examples are of two dimensional form, the skilled person would be readily aware that the present invention is applicable to three dimensional approaches.
The generated model consists of values of the coefficient Vp and, possibly, other physical values or coefficients, typically defined over a discrete grid representing the body region to be imaged.
The accuracy of the starting model is of significant importance to successfully image regions of the body that include bone and/or gas. It is of particular significance with regard to imaging of the head since the skull is a closed system and the measurements depend heavily upon the properties of the skull. Therefore, a starting model with an effective starting configuration for the skull is critical for successful convergence.
Therefore, it is necessary to ensure that the starting model is sufficiently representative of the region of the skull through which imaging is to take place in order to ensure successful imaging of the brain and soft tissue within.
The step of providing a starting model that is representative of the skull involves the selection, generation and/or modification of a starting model for the skull for use with the imaging method. These alternatives will now be described.
The skull component of the starting model may be generated in a number of different, non-limiting, ways. However, they can broadly be placed in two categories. The first is where one or more starting model skull components are generated in advance of the imaging procedure and form a database of starting model skull components. An appropriate skull component is then selected based on an empirical parameter relating to the subject or a measurement carried out on the subject.
The second is that the skull component of the starting model is procedurally generated or modified in response to an empirical parameter or measurement of the subject.
Firstly, it is necessary to acquire data in advance of the imaging measurement in order to generate one or more starting models. Key to the process is an accurate starting model of the skull.
This may be done in numerous ways. For example, reflection ultrasound methods could be used to obtain a starting model for the skull. This may involve, for example, using high frequencies or focused beams.
Alternatively or additionally, low intensity X-ray computed tomography (CT) scans could be used to provide an estimation of the skull properties which can be converted to velocities to generate a starting model. This data can be acquired beforehand either on the subject (e.g. in a screening process) or from a range of different subjects. This method can be done for bone or in the presence of air.
Alternatively, data may be obtained from other sources (e.g. MRI). As a further alternative, low-frequency transmission and reflection ultrasound could be used with FWI to obtain a full model of the skull. The low frequencies mean that cycle skipping can be avoided to build a better skull component of the starting model which is closer to the global minimum.
However, the method of acquisition of the starting model skull component is not material to the present invention. What is of benefit is that one or more skull components of the starting model are available to be selected from or modified in response to empirical data or measurements of the subject to be imaged.
The next stage is to acquire data relating to the subject to be imaged. This may be done in any suitable manner. The data acquired may relate to physical measurements made on the subject in situ, general characterising empirical data such as age, sex, weight, height or a combination of the above.
If measurements are made on the subject, this is to determine the structure and morphology of at least the skull. Key data to be acquired may relate to the thickness, composition and morphology.
The measurements made on the subject may include the observed data set acquired in step 100 as will be described below. Alternatively or additionally, other measurements may be used. For example, shear sensors could be placed on the skull and the skull excited with a shear stress. The surface waves propagating along the skull could then be measured and used to obtain information about the thickness and morphology of the skull.
Alternatively or additionally, lasers and/or direct physical measurements could be used to measure the outer geometry of the head.
As a further alternative or addition, the observed data set acquired in step 100 could be used to inform the choice of starting model. For example, in the case that a database of predetermined starting models has been provided, a corresponding predicted data set for each starting model could also be provided. Then, at least a part of the observed data set acquired in step 100 could be matched with the starting model predicted data sets to find the closest match. The starting model which represents the closest match could then be used in the full imaging process.
As a further alternative or addition, a part of the observed data set, for example at early travel times (representing reflections or refractions from the skull), could be used to generate or modify a model of the skull to generate a more accurate starting model for the full inversion process.
Any of the above empirical data or measured data parameters could be used to select from a pre-determined group of starting model skull components or to modify or generate a starting model skull component which has suitable accuracy to enable the imaging method to be carried out.
Once the starting model skull component has been generated, then a soft tissue component can be added to form the full starting model for the method of the present invention.
Since, as described above, the accuracy of the soft tissue component of the starting model is less critical than that of the skull component, the soft tissue component may comprise simply a homogeneous layer having an acoustic velocity around 1500 m/s (which is typical of soft tissue). Alternatively, other configurations of soft tissue component may be provided to account for the acoustic velocity of soft tissue varying from around 1450 m/s for fat to about 1730 m/s for skin. However, it is generally not necessary to provide for the variations in morphology and structure as is required from the skull component.
In summary, then, the starting model comprises a first component including elements having an acoustic velocity in excess of 2300 m/s (in the case of a skull), and a second component comprising elements having an acoustic velocity within the range of 1400-1750 m/s. This is the case for a starting model of the head as described in this non-limiting embodiment. However, other values may be used for the two components if other body parts are to be imaged (described later)
The method then proceeds to step 104.
Step 104: Generate predicted data set
At step 104, a predicted ultrasonic data set is generated. The predicted data is required to correspond to the same source-receiver location data positions as the actual measured trace data from the ultrasonic measurement so that the modelled and observed data can be compared. In other words, the predicted data set corresponds discrete point to discrete point to the observed data set. The predicted data set is generated for the same measurement parameter(s) at the same frequency or frequencies.
Predicted ultrasonic data may be generated based on an analysis of the acoustic isotropic two-way wave equation as set out below in equation 17):
where the acoustic pressure p and driving source s vary in both space and time, and the acoustic velocity c and density p vary in space. This equation applies to small-amplitude pressure waves propagating within an inhomogeneous, isotropic, non-attenuating, non-dispersive, stationary, fluid medium. It is relatively straightforward to add elastic effects, attenuation and anisotropy to the wave equation. Introduction of these parameters changes the detailed equations and numerical complexity, but not the general approach.
The wave equation 17) represents a linear relationship between a wave field p and the source s that generates the wave field. After discretisation (with, for example, finite differences) we can therefore write equation 17) as a matrix equation 18):
Ap=s 18)
where p and s are column vectors that represent the source and wavefield at discrete points in space and time, and A is a matrix that represents the discrete numerical implementation of the operator set out in equation 3):
Although the wave equation represents a linear relationship between p and s, it also represents a non-linear relationship between a model m and wavefield p. Thus equation 17) can be rewritten as equation 20):
G(m)=p 20)
where m is a column vector that contains the model parameters. Commonly these will be the values of c (and p if density is an independent parameter) at every point in the model, but they may be any set of parameters that is sufficient to describe the model, for example slowness 1/c, acoustic modulus c2p, or impedance cp.
In equation 20), G is not a matrix. Instead it is a non-linear Green's function that describes how to calculate a wavefield p given a model m.
From the above analysis, predicted ultrasonic data can be generated for one or more physical parameters in the time domain. If done in the frequency domain it could be done for one or more selected frequencies. This forms the predicted ultrasonic data set dpred. The method now proceeds to step 106.
Step 106: Construct misfit function
At step 106, a misfit function is configured. In one example, the misfit function (or objective function) is configured to measure the dis-similarity between the observed and predicted data sets. Alternatively, an objective function may be configured that measures similarity; in this case, step 108 will operable to maximise rather than minimise the objective function.
The method proceeds to step 108.
Step 108: Minimise or maximise the misfit function
This may be done by any suitable method. In this embodiment, the gradient method is used to minimise the misfit between the two data sets.
The method then proceeds to step 110.
Step 110: Update model
At step 110, the model is updated using the gradient obtained in step 108. The model update derives from the gradient of equation 11) i.e. the partial derivative with respect to a point perturbation of the model m at each position. Ultimately gradients from separate sources will summed when forming the final model update.
As with the computational structure of conventional FWI method as used in seismic modelling, this is the product of two wavefields: an incident wavefield emitted by a source at a source location and a back-propagated wavefield which is emitted by a (multi-point) source located at the receiver positions.
As noted above, gradient methods can be enhanced using an approximate form for the Hessian and conjugate directions for the model update. Further a step-length is then calculated to scale the search direction and give the final model update.
For any useful measure of misfit or similarity the predicted ultrasonic data set will move towards the observed ultrasonic data set. Thus, the starting model will move towards the true model. The method now proceeds to step 112.
Step 112: Convergence criteria met?
At step 112 it is determined whether convergence criteria have been met. For example, the method may be deemed to have reached convergence when the difference between the data sets reaches a threshold percentage or other value. If the criteria as set out above have been met, then the method proceeds to step 114 and finishes with the resultant model generated. If the criteria have not been met, then the method proceeds back to repeat steps 104 to 110 as discussed above.
When, at step 114, it has been determined that the convergence criteria has been met, the method finishes and the modelled region of the subject's body is deemed to be sufficiently accurate to be used for medical imaging.
Information can then be extracted from the imaged region of the subject's body. For example, the imaged region may be used to identify a pathology or to assist in the diagnosis of a pathology.
If applied to the head, the imaging method may be used to detect the presence and/or to determine the absence of congenital and/or acquired abnormalities in tissue composition and/or morphology within a region or regions within the intra-cranial cavity.
The images may also be used, for example to detect, identify, characterise, locate and/or differentiate between normal tissue and abnormal tissue caused by intra-cranial pathologies of vascular supply to the brain involving an interruption to, and/or a reduction of, blood flow to brain tissue and/or extradural, subdural, subarachnoid, intra-cerebral and/or intra-ventricular haemorrhage.
Alternatively, if the imaged region of the body is one containing bones then the imaging method could be used to detect, identify, characterise and/or locate any abnormality due to muscular or osseus injuries.
If the torso is the subject of the imaging, the method could be used to identify the totality of the viscera to detect, identify, characterise and/or locate any abnormality, for instance, due to a tumour in the pancreas or kidney stones.
If a pregnant woman is imaged, the recovered acoustic properties of the abdominal part of a pregnant woman to detect, identify, characterise, locate and/or image a foetus.
Alternatively, the method can be used to detect, identify, characterise and locate tumours that are otherwise obscured by the presence of bone and/or gas within the region being imaged.
Whilst the embodiment of the method above has been described with specific references to the head and skull, the method is applicable to other regions of the body which may contain interfaces between bone, gas and/or soft tissue and which cannot be imaged by conventional means.
The described method of the present invention can be generalised to other parts of the body where a starting model is built up and comprises at least two components—a first component comprising a plurality of model parameters representative of the physical properties and morphology of bone and/or gas within the body part of the subject to be imaged and having at least one modelled region having an acoustic velocity below 700 m/s and/or above 2300 m/s, and a second component comprising a plurality of parameters representative of the physical properties of the soft tissue within the body part of the subject to be imaged. In general, the soft tissue component of the starting model comprises elements having an acoustic velocity within the range of 1400-1750 m/s.
In some cases, a variation to the described embodiment of the method may be required. For example, the method of obtaining a starting model may be different for certain regions of the body when compared to the head. It may, for instance, not be possible to maintain a database of different starting model components in cases where there may be large physical, genetic and/or geometric differences between subjects.
For example, the imaging of a joint may be problematic since it is challenging to provide a range of starting models which capture the possible variations in joint angle, size, construction etc. However, imaging of such regions is generally not as time-sensitive as imaging of the head and brain (e.g. in the case of stroke). Therefore, a measurement of the structure of the body region may be made simultaneously or sequentially using, for example, a low intensity X-ray CT image from which a starting model is subsequently derived.
The skilled person would be readily aware of the suitable environments in which data could be gathered for imaging and analysis purposes as set out in the present disclosure.
In aspects, the embodiments described herein relate to a method of extracting information from a digital image. However, the embodiments described herein are equally applicable as an instruction set for a computer for carrying out said method or as a suitably programmed computer.
The methods described herein are, in use, executed on a suitable computer system or device running one or more computer programs formed in software and/or hardware and operable to execute the above method. A suitable computer system will generally comprise hardware and an operating system.
The term ‘computer program’ is taken to mean any of (but not necessarily limited to) an application program, middleware, an operating system, firmware or device drivers or any other medium supporting executable program code.
The term ‘hardware’ may be taken to mean any one or more of the collection of physical elements that constitutes a computer system/device such as, but not limited to, a processor, memory device, communication ports, input/output devices. The term ‘firmware’ may be taken to mean any persistent memory and the program code/data stored within it, such as but not limited to, an embedded system. The term ‘operating system’ may taken to mean the one or more pieces, often a collection, of software that manages computer hardware and provides common services for computer programs.
The comparison step may also be conducted making use of previous measurements on a healthy or diseased population of reference joints for which values or average values are stored in a database or memory location in such a computer. The computer may be programmed to display the results of the comparison as a read out.
The methods described herein may be embodied in one or more pieces of software and/or hardware. The software is preferably held or otherwise encoded upon a memory device such as, but not limited to, any one or more of, a hard disk drive, RAM, ROM, solid state memory or other suitable memory device or component configured to software. The methods may be realised by executing/running the software. Additionally or alternatively, the methods may be hardware encoded.
The method encoded in software or hardware is preferably executed using one or more processors. The memory and/or hardware and/or processors are preferably comprised as, at least part of, one or more servers and/or other suitable computing systems.
Number | Date | Country | Kind |
---|---|---|---|
1621436.3 | Dec 2016 | GB | national |
Number | Date | Country | |
---|---|---|---|
Parent | 16468029 | Jun 2019 | US |
Child | 18374037 | US |