IMAGE RECONSTRUCTION METHOD BASED ON A TRAINED NON-LINEAR MAPPING

Abstract

The present invention concerns an array processing image reconstruction method comprising: receiving a set of echo waveforms as measurements from the object of interest; defining a measurement model linking an unknown image of the object to the measurements; defining a data fidelity functional using the measurement model; defining a regularisation functional using a trained non-linear mapping, the regularisation functional comprising prior knowledge about the unknown image; defining an optimisation problem involving the data fidelity functional and the regularisation functional for obtaining a first image estimate of the unknown image; and solving the optimisation problem to obtain the first image estimate.

Claims

1. An image reconstruction method comprising: receiving a set of waveforms as measurements y from an object of interest; defining a measurement model H linking an unknown image x of the object to the measurements y; determining a data fidelity functional using the measurement model, the data fidelity functional measuring a distance between H(x) and the measurements y; determining a regularisation functional using a trained non-linear mapping f.sub.θ, the mapping being non-linear with respect to its input data, the regularisation functional comprising prior knowledge about the unknown image; defining an optimisation problem involving the data fidelity functional and the regularisation functional for obtaining a first image estimate {circumflex over (x)} of the unknown image; and solving the optimisation problem to obtain the first image estimate {circumflex over (x)}.

2. The method according to claim 1, wherein the method comprises, prior to defining the data fidelity functional and the regularisation functional , training the non-linear mapping by using a training data set.

3. The method according to claim 2, wherein the training data set comprises images acquired experimentally or by numerical simulation.

4. The method according to claim 1, wherein the trained non-linear mapping is a trained artificial neural network.

5. The method according to claim 1, wherein the optimisation problem is expressed as $\begin{array}{cc}\hat{x}=\underset{x\in {ℝ}^N}{\mathrm{argmin}}\left\{\mathrm{��}\left(x\right)+\mathrm{\lambda ℛ}\left(x\right)\right\},& \end{array}$ where λ∈.sup.+ is a regularisation parameter controlling the weighting of the regularisation functional .

6. The method according to claim 1, wherein the data fidelity functional and/or the regularisation functional has at least one of the following properties: convexity, non-convexity, differentiability and non-differentiability.

7. The method according to claim 1, wherein solving the optimisation problem comprises executing an iterative algorithm.

8. The method according to claim 7, wherein the iterative algorithm used for solving the optimisation problem is a gradient descent method or a generic projection-based method.

9. The method according to claim 7, wherein the method comprises setting an iteration step size to one, setting an initialisation parameter to zero and including only one iteration in the iterative algorithm and then applying the non-linear mapping directly to a second image estimate {tilde over (x)} for obtaining the first image estimate {circumflex over (x)} of the object of interest, wherein the second image estimate {tilde over (x)} has a lower quality in terms of image artefacts than the first image estimate {tilde over (x)}.

10. The method according to claim 7, wherein the trained non-linear mapping varies from one iteration step to another.

11. The method according to claim 7, wherein the method further comprises transmitting a set of waveforms to an object of interest.

12. The method according to claim 11, wherein the number of transmitted waveforms is between one and five, and more specifically between one and three, and more specifically one.

13. The method according to claim 1, wherein the image reconstruction method is an ultrasound image reconstruction method.

14. The method according to claim 1, wherein the regularisation functional is defined as one of the following: (x)=½x.sup.T [x−f.sub.θ(x)]; R(x)=½∥x−f.sub.θ(x)∥.sub.2.sup.2; (x)=½∥r.sub.θ(x)∥.sub.2.sup.2, where r.sub.θ: .sup.N.fwdarw..sup.N is a non-linear mapping trained to predict the negative noise to be applied to x; and (x)=½∥e.sub.θ.sub.e(x)∥.sub.2.sup.2, where e.sub.θ.sub.e is an encoder with trainable parameters θ.sub.e.

15. An image reconstruction apparatus being configured to perform operations comprising: receive a set of waveforms as measurements y from an object of interest; define a measurement model H linking an unknown image x of the object to the measurements; determine a data fidelity functional using the measurement model, the data-fidelity functional measuring a distance between H(x) and y; determine a regularisation functional using a trained non-linear mapping, the mapping being non-linear with respect to its input data, the regularisation functional comprising prior knowledge about the unknown image; define an optimisation problem involving the data fidelity functional and the regularisation functional for obtaining a first image estimate {circumflex over (x)} of the unknown image; and solve the optimisation problem to obtain the first image estimate {circumflex over (x)}.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Other features and advantages of the invention will become apparent from the following description of a non-limiting example embodiment, with reference to the appended drawings, in which:

[0015] FIG. 1 is a block diagram illustrating some elements of an imaging system, which are useful for understanding the teachings of the present invention, according to an example of the present invention;

[0016] FIG. 2 is a flow chart illustrating the proposed image reconstruction method according to an example of the present invention;

[0017] FIG. 3 is a block diagram illustrating a direct non-linear mapping according to an example of the present invention;

[0018] FIG. 4 is a block diagram illustrating a residual non-linear mapping according to an example of the present invention;

[0019] FIG. 5 is a block diagram illustrating a residual non-linear mapping having an “encoder-decoder” architecture according to an example of the present invention;

[0020] FIG. 6 is a flow chart illustrating the proposed image reconstruction method using a gradient descent algorithm according to an example of the present invention;

[0021] FIG. 7 is a flow chart illustrating the proposed image reconstruction method using a generic projection-based algorithm according to an example of the present invention; and

[0022] FIG. 8 is a block diagram illustrating an example deep residual convolutional neural network, which may be used in the present invention.

DETAILED DESCRIPTION OF AN EMBODIMENT OF THE INVENTION

[0023] An embodiment of the present invention will now be described in detail with reference to the attached figures. This embodiment is described in the context of pulse-echo US imaging, but the teachings of the invention are not limited to this environment. Identical or corresponding functional and structural elements which appear in the different drawings are assigned the same reference numerals. In the present description, vector variables and functions outputting a vector are denoted with a bold font.

[0024] FIG. 1 is a simplified block diagram showing some elements of the imaging system or apparatus 1, which are useful for understanding the teachings of the present invention. The proposed imaging system, as shown in FIG. 1 comprises a US probe 3 capable by means of a transmit and receive module 3 to transmit an acoustic beam in a medium of interest and receive backscattered echoes. The images are reconstructed using the image reconstruction module 7 which takes advantage of a trained non-linear mapping or trained non-linear function used in conjunction with an image reconstruction algorithm or process carried out by the image reconstruction module 7. The trained non-linear mapping is non-linear with respect to its input data (in this example images) and is arranged to map an input vector space to an output vector space. A post-processing and display module 9 is arranged to post-process the reconstructed image. The main advantage of the proposed image reconstruction method is its capability of reconstructing high-quality images, free of common array processing artefacts, with a number of transmit-receive events dramatically reduced compared to the number of transmit-receive events required to obtain a similar image quality with state-of-the-art image reconstruction methods.

[0025] The US probe 3 in this example has a finite aperture comprising a linear array of piezoelectric elements, also referred to as transducers. Such a probe is suitable for the purpose of imaging in 2D, and may include a cylindrical focusing scheme to concentrate the acoustic energy, in elevation, in the vicinity of the image plane. However, this kind of probe is merely an example. The proposed method could be easily adapted to different probe geometries (such as convex probes, phased arrays etc.) or technologies (including polyvinylidene fluoride (PVDF) copolymers and capacitive micromachined ultrasonic transducers (CMUT) etc.). Similarly, the proposed method can be extended to cases using 2D matrix probes, designed to transmit acoustic beams and collect echoes from the insonified volume in order to provide volumetric imaging, whereby the imaging process produces a 3D representation of the volume under scrutiny. In 2D imaging configurations, probes with transducer elements arranged along a straight planar, convex or concave line are conventionally used. In 3D imaging configurations, matrix-array probes with transducer elements arranged on a planar, convex or concave surface are conventionally used. These transducers elements are capable of converting electrical excitation into acoustic pressure waveforms and record acoustic pressure waveforms into electrical signals. This property is used for both transmitting acoustic beams in the medium of interest and receiving the corresponding backscattered echoes, which arise from local acoustic impedance variations. The probe 3 thus forms an array processing configuration having an array structure comprising a set of transducers, which can act as sensors or actuators, that are spatially separated.

[0026] The transmit and receive module 5 in this example comprises a transmit path and a receive path and is connected to the US probe 3. The transmit path generally contains frontend electronics, which may include an electrical excitation circuitry with controlled timing, voltage-amplitude and waveform shaping means. It is used for controlling and shaping the desired transmit beams. The receive path generally contains electrical circuitry with associated programmable gain controls (such as fixed gain and time-gain compensation) for counteracting the object (e.g. tissue) attenuation, and analogue-to-digital conversion means followed by digital signal processing means. The received signals or measurements are then sent to the image reconstruction module 7.

[0027] The objective of the image reconstruction module 7 is to reconstruct an image of the medium under scrutiny from the received signals. To this end, the image reconstruction module 7 implements an image reconstruction method, process or algorithm, which uses a trained non-linear mapping defined as f.sub.θ. The trained non-linear mapping may be an ANN whose trainable parameters θ have been optimised to reduce or even remove common US imaging artefacts. In other words, the non-linear mapping may be a machine learning configuration having trainable parameters. The training of these trainable parameters is carried out by using optimisation algorithms using a dedicated dataset. Once trained, the trained non-linear mapping is used in the image reconstruction method to produce high-quality images with a minimal number of transmit-receive events.

[0028] The image reconstruction method solves an ill-posed inverse problem of the form

y=H(x)+n,  (1)

where x∈.sup.N ( denoting real numbers) is the unknown image with N being the number of image elements, i.e. pixels or voxels composing the image, y∈.sup.M are the measurements with M being the number of measurement samples, H: .sup.N.fwdarw..sup.M is the measurement model, which accounts for the physical measurement process (which in the case of US may take advantage of a relation between time-of-arrival and distance travelled by the waveform, as determined by the speed-of-sound in the medium), and n∈.sup.M is the measurement noise. Finding the unknown image x or its estimate {circumflex over (x)} can be performed using an iterative algorithm to solve an optimisation problem of the form:

[00001]$\begin{array}{cc}\hat{x}=\underset{x\in R^N}{\mathrm{argmin}}\left\{\mathrm{��}\left(x\right)+\mathrm{\lambda ℛ}\left(x\right)\right\},& \left(2\right)\end{array}$

where : .sup.N.fwdarw..sup.+ is the data fidelity functional or term (also known as a data discrepancy functional), : .sup.N.fwdarw..sup.+ (R.sup.+ denoting positive real numbers) is the regularisation functional or term commonly referred to as a prior encoding prior knowledge about the unknown image x, and λ∈.sup.+ is a hyper-parameter controlling the weighting of the regularisation functional with respect to the data fidelity functional. A functional is understood to refer to a (e.g. linear) mapping from a vector space into the real numbers . Instead of using standard image processing priors, such as sparsity in fixed bases or frames, which are of limited appeal in US imaging, a highly efficient and representative data-driven regularisation functional can be constructed using a trained non-linear mapping f.sub.θ. Many well-known iterative algorithms are available in the literature to solve the problem of Equation (2), such as gradient descent, proximal gradient descent and projected gradient descent, depending on the mathematical properties of both functionals and . Such algorithms may benefit from the differentiability of the regularisation functional , which implies the non-linear mapping f.sub.θ to be differentiable. In such a setting, the use of a trained NN as the non-linear mapping f.sub.θ may be beneficial since they can be built differentiable or can exploit standard techniques to accurately approximate their derivative. Moreover, efficient evaluation of their derivative can be readily exploited to speed up the computations involved in each iteration of the algorithm.

[0029] The flow chart of FIG. 2 summarises the proposed image reconstruction method according to an example. In step 11, the probe 3 transmits by using a set of transmitting elements a set of waveforms in a medium to be reconstructed at least partially. However, it is to be noted that step 11 is optional since the object of interest could simply be exposed to a physical excitation, which may include US (acoustic waveforms), thermal, or other excitations, such as electromagnetic-radiation pulses in photoacoustic imaging. In step 13, the probe 3 or more specifically its sensor(s) receives or acquires by using a set of receiving elements a set of (echo) waveforms arising from the physical scattering of the insonified medium. In other words, in this step a set of waveforms (or signals) is received or acquired as measurements from the object of interest stimulated by (or exposed to) a physical excitation. In step 15, a measurement model is defined, determined or identified linking the unknown image of the medium to the recorded set of measurements. In step 17, a data fidelity functional is constructed, determined or identified using the defined measurement model as will be explained later. In step 19, a data-driven regularisation functional or prior is constructed, determined or identified using a trained non-linear mapping f.sub.θ as will be explained later. In step 21, an optimisation problem involving the data fidelity functional and the regularisation functional is defined, determined or identified for obtaining a meaningful or high-quality estimate {circumflex over (x)} of the unknown image. In step 23, the optimisation problem is solved for instance by means of an iterative algorithm or solver until a convergence criterion is reached to obtain the image estimate {circumflex over (x)} of the medium under scrutiny. The algorithm may comprise between 10 and 100 iterations or more specifically between 10 and 50 iterations or between 1 and 20 iterations. The steps 15 to 21 may be carried out by the image reconstruction module 7.

[0030] After image reconstruction, the reconstructed data (i.e. an image) are sent to a post-processing and display module 9 as further shown in FIG. 1. The post-processing step covers a range of applications, such as: [0031] B-mode imaging; [0032] scan conversion; [0033] colour Doppler imaging; [0034] vector flow imaging; and [0035] elastography.

[0036] For B-mode imaging using RF data, envelope detection is applied onto the reconstructed data (i.e. the image). The envelope detection can, for instance, be achieved by means of the Hilbert transform followed by a magnitude detection and optional low-pass filtering. It can also be achieved by squaring and low-pass filtering the signal. In the case of B-mode imaging using IQ data, the magnitude of the complex signal is extracted. The envelope detection step is followed by a gain adjustment and dynamics compression steps. For Doppler and elastography, the reconstructed RF or IQ data (i.e. the image) are used directly, without post-processing.

[0037] In the above description, a brief overview of the invention was given. The invention will be next described in more detail with reference to some concrete examples.

[0038] US imaging, and more generally array processing, suffers from artefacts due to several factors, some of which are described below. Typical US probes have a finite aperture which is subdivided into an array of transducer elements. Each of these transducer elements has bandlimited, or bandpass, electro-acoustic and acousto-electric impulse responses allowing them to convert electrical excitation into acoustic pressure and record acoustic pressure into an electrical signal. This implies that an ideal point reflector in tissue can only be recorded as a bandlimited, or bandpass, RF pulse waveform. The finite dimension nature of transducer elements makes them diffracting objects with directional sensitivity. The probe aperture, the arrangement of these transducer elements within the aperture, as well as the delays and amplitudes applied for shaping the transmitted beam are governing factors for the array transmit and receive spatial impulse responses (SIRs). The number and type of transmitted beams and corresponding DAS beamforming operations are key factors in the resulting image quality. The probe, the arrangement of transducer elements, their shapes and physical properties together with the transmit and receive settings are defining parameters of a US imaging configuration. Such an imaging configuration is coupled to an image reconstruction method, such as the well-known DAS, to produce an image from the recorded measurements. A standard quality measure of a resulting image is given by the point-spread-function (PSF). Wave propagation effects such as scattering and attenuation also impact the overall image quality. The entire physical imaging process entails the appearance of imaging artefacts, caused most notably by side-lobe and grating-lobe levels, attenuation, reverberation and multiple scattering, which are visualised as reduced contrast, ghost images in incorrect locations, non-physical striped image structures, shadowing, etc.

[0039] As an example, let us define two imaging configurations as detailed in Table 1. A low-quality configuration, called LQ, uses a linear probe composed of 128 transducer elements and is centred at 5.133 MHz with 67% bandwidth. Images are reconstructed from a single PW measurement, i.e. all the elements of the array transmit an acoustic wave in the medium at the same time, resulting in an unfocused quasi-plane wavefront. The backscattered echoes are measured by all the elements of the array from which the image is reconstructed by DAS at every pixel location of the image estimate. A high-quality configuration, called HQ, uses a similar probe with twice as many transducer elements within the same aperture. Images are reconstructed from 256 SA measurements, which can be considered as a gold standard in US imaging with a given aperture. SA measurements are obtained as follows. A single element of the array transmits an unfocused wavefront in the medium. The backscattered echoes are measured by all the elements of the array from which pixels of an intermediate RF image are reconstructed by DAS. This procedure is repeated for each element of the array and the final RF image is obtained by averaging all intermediate RF images.

TABLE-US-00001 TABLE 1 Parameter LQ HQ Element number 128 256 Pitch 300 μm 150 μm Element width 270 μm 150 μm Aperture 38.1 mm 38.1 mm Element height 5 mm 5 mm Elevation focus 20 mm 20 mm Centre frequency 5.133 MHz 5.133 MHz Bandwidth 67% 67% Transmit frequency 5.208 MHz 5.208 MHz Excitation cycles 2.5 2.5 Sampling frequency 20.832 MHz 20.832 MHz Transmit-receive scheme 1 PW 256 SA

[0040] By assessing the image reconstruction results, it can be noted that the PSF of the LQ configuration suffers far more from artefacts than the HQ configuration. The side lobes and ghost artefacts, which mostly impact the image contrast, are drastically reduced in the HQ configuration. The main lobe, which has a direct impact on the image resolution, is tightened in the HQ configuration. The grating lobes, which can have devastating blurring and shadowing effects, are present in the LQ configuration. This results from the fact that the considered array-pitch is greater than half the wavelength. These grating lobes are entirely removed in the HQ configuration, which has a half-wavelength pitch. Overall the reconstructed images when using the HQ configuration have a much higher quality than the images reconstructed from the LQ configuration.

[0041] Let W⊂.sup.N be the subspace of low-quality images reconstructed from an imaging configuration inducing a substantial number of artefacts, where N represents the dimension of an image, such as the number of pixels in two-dimensional (2D) space or the number of voxels in three-dimensional (3D) space. Let V⊂.sup.N be the subspace of high-quality images. Such high-quality images may be reconstructed by an optimal US imaging configuration which reduces or even entirely removes these artefacts. It may also be high-quality images acquired with another imaging modality immune to these artefacts or the exact physical medium properties if available. We define the mapping F: V.fwdarw.W which transforms a high-quality image or a high-quality image estimate (also referred to as a first image estimate) {circumflex over (x)}∈V to a corresponding low-quality image or a low-quality image estimate (also referred to as a second image estimate) {tilde over (x)}∈W such that

F({circumflex over (x)})={tilde over (x)}.  (3)

[0042] Many different mappings F may exist, depending the subspaces V and W, but all of them are ill-posed in essence. Indeed, since the level of artefacts present in low-quality images is higher than the level of artefacts present in high-quality images, it implies that F leads to a loss of information. In other words, the mapping F acts as a blurring operator onto the high-quality image {circumflex over (x)}. Hence, finding the inverse mapping to recover a high-quality image {circumflex over (x)} from its corresponding low-quality image {tilde over (x)} is not possible by means of state-of-the-art methods.

[0043] As an example, let us consider the HQ and LQ configurations defined in Table 1. In order to formalise the resulting image subspaces W and V, let us define Y.sub.L⊂.sup.M.sup.L and Y.sub.H⊂.sup.M.sup.H as the subspaces of measurements acquired with the LQ and HQ configurations, where M.sub.L and M.sub.H are their corresponding number of measurement samples. In this example, it can be noted that M.sub.H is substantially greater than M.sub.L since the HQ configuration is composed of twice the number of transducer elements and requires 256 times more transmit-receive events than the LQ configuration. Equipped with these definitions, one can express a low-quality image {tilde over (x)} and a high-quality image {circumflex over (x)} as,

{tilde over (x)}=D.sub.L(y.sub.L),y.sub.L∈Y.sub.L,  (4)

{circumflex over (x)}=D.sub.H(y.sub.H),y.sub.H∈Y.sub.H,  (5)

where D.sub.L:Y.sub.L.fwdarw.W and D.sub.H:Y.sub.H.fwdarw.V are the imaging operators for the LQ and HQ configurations, respectively. Hence, the subspace of low-quality images can be expressed as W=D.sub.L(Y.sub.L) and the subspace of high-quality images as V=D.sub.H(Y.sub.H). In the particular example described above, DAS is used as the image reconstruction method, hence both operators D.sub.L and D.sub.H are linear operators. It is interesting to note that the reconstruction method used to reconstruct an image from the set of SA measurements obtained with the HQ configuration is in fact a process which averages coherently multiple views of the same image. Since the images reconstructed with the HQ configuration are obtained with far more views than images reconstructed from the measurements obtained with the LQ configuration, they inherently contain more information. Therefore, it is possible to express the mapping F as defined in Equation (3), which, in this example, is linear, e.g. by means of principal component analysis (PCA). However, it is clear that finding its inverse mapping by means of state-of-the-art methods is not possible since it amounts to recovering lost information.

[0044] In order to circumvent this limitation, it is common to formulate an imaging reconstruction task as an inverse problem. Let us define an unknown image x∈.sup.N and the measurements y∈.sup.M sensed by a US imaging configuration, where M represents the number of measurement samples recorded by the transducer elements of the probe. The corresponding inverse problem can be written as

y=H(x)+n,  (6)

where H: .sup.N.fwdarw..sup.M is the measurement model or forward model (having x as one of its inputs) accounting for the physical measurement process and n∈.sup.M is the measurement noise. It is important to note that the measurement model may be non-linear. Even though the proposed formulation is expressed with a measurement model taking values in .sup.N (real numbers), it can be readily extended to the case where the measurement model takes values in .sup.N (complex numbers). In the special case of a linear inverse problem the measurement model can be expressed as a matrix H∈.sup.M×N such that H(x)=Hx. In this context, the problem of Equation (6) becomes

y=Hx+n.  (7)

[0045] Since the problem of Equation (6) is generally ill-posed, it may be solved by formulating the unconstrained optimisation problem

[00002]$\begin{array}{cc}\hat{x}=\underset{x\in {ℝ}^N}{\mathrm{argmin}}\left\{\mathrm{��}\left(x\right)+\mathrm{\lambda ℛ}\left(x\right)\right\}.& \left(8\right)\end{array}$

[0046] In the case where both functionals and are convex, the problem of Equation (8) can be equivalently formulated as the constrained minimisation problem

[00003]$\begin{array}{cc}\hat{x}=\underset{x\in C_{ℛ}}{\mathrm{argmin}}\mathrm{��}\left(x\right),& \left(9\right)\end{array}$

where ={x∈.sup.N|(x)≤τ} is a closed convex set and T is uniquely defined by the regularisation parameter λ.

[0047] The data-fidelity functional , which measures the distance between H(x) and the measurements y, is often expressed as the squared .sub.2-norm

(x)=½∥H(x)−y∥.sub.2.sup.2  (10)

[0048] Different regularisation functionals , which are often called priors, are commonly used in the literature depending on the nature of the unknown image x. Some popular priors are the squared .sub.2-norm, .sub.1-norm, wavelet-sparsity or total variation (TV). Even though such priors have demonstrated their efficacy in standard image processing tasks, they remain far from optimal in the context of US imaging due to the highly specific nature of US images. Learned dictionaries combined with a suitable norm are also commonly used as regularisation functionals. However, they are computationally demanding and only applicable to patches of an image, clearly limiting the appeal of such techniques.

[0049] The present invention proposes a strategy which involves using a suitably trained non-linear mapping f.sub.θ:.sup.N.fwdarw..sup.N with trainable parameters θ to construct a highly effective data-driven prior for the unknown image x. FIG. 3 shows a simplified block diagram of such a non-linear mapping f.sub.θ. The trained non-linear mapping can typically be a deep neural network (DNN) whose trainable parameters θ have been trained to reduce or even remove common US imaging artefacts.

[0050] Equipped with the trained non-linear mapping f.sub.θ, a data-driven prior can be expressed through a regularisation functional . A typical example of may be defined as the inner product between the image x and the residual term [x−f.sub.θ(x)] such that

(x)=½x.sup.T[x−f.sub.θ(x)],  (11)

where x.sup.T represents the transpose of x. In the case of a complex vector x∈.sup.N, the transpose operation should be replaced by the conjugate transpose (or Hermitian transpose) denoted x.sup.H.

[0051] Another way of defining such a regularisation functional may be the squared .sub.2-norm of said residual, which is explicitly expressed as

(x)=½∥x−f.sub.θ(x)∥.sub.2.sup.2.  (12)

[0052] It is clear that in both examples given in Equations (11) and (12), (x) will tend to zero as the image x tends to f.sub.θ(x). In some scenarios, especially when considering DNNs, it may be desirable to design a non-linear mapping f.sub.θ as a residual mapping such that

f.sub.θ(x)=x+r.sub.θ(x),  (13)

where r.sub.θ: .sup.N.fwdarw..sup.N is a non-linear mapping trained to predict the negative noise to be applied to x.

[0053] FIG. 4 shows a block diagram of such a residual non-linear mapping. In this case, the regularisation functional defined in Equation (12) can be simplified for example to become

(x)=½∥r.sub.θ(x)∥.sub.2.sup.2.  (14)

[0054] It is clear that, in this case, (x) will tend to zero as the predicted negative noise r.sub.θ(x) tends to zero.

[0055] Another popular NN architecture is the “encoder-decoder” architecture, where an encoder e.sub.θ.sub.e with trainable parameters θ.sub.e first maps the input image x to a latent space such that e.sub.θ.sub.e: .sup.N.fwdarw..sup.E, where E represents the dimension of the latent space, which is usually smaller than the dimension of the image space N. The decoder d.sub.θ.sub.d with trainable parameters θ.sub.d then maps the latent space back to the image space such that d.sub.θ.sub.d: .sup.E.fwdarw..sup.N. In its residual form, such an architecture can be expressed as

f.sub.θ(x)=x+r.sub.θ(x)=x+d.sub.θ.sub.d(e.sub.θ.sub.e(x)),  (15)

where the trainable parameters θ include both θ.sub.e and θ.sub.d. FIG. 5 shows a block diagram of a residual non-linear mapping with an “encoder-decoder” architecture. For such a trained non-linear mapping f.sub.θ, one may express the regularisation functional as

R(x)=½∥e.sub.θ.sub.e(x)∥.sub.2d.sup.2,  (16)

which might be advantageous in terms of computational complexity.

[0056] Equipped with the data-driven regularisation functional expressed in terms of a non-linear mapping f.sub.θ appropriately trained for a given application, one can use it in a suitable image reconstruction method to produce high-quality images from a US imaging configuration only capable of producing low-quality images with a traditional image reconstruction method such as the well-known DAS. The proposed image reconstruction method typically requires solving and optimisation problem such as Problem (8) by means of an iterative algorithm. Many algorithms may be deployed for solving Problem (8) depending on the mathematical properties, such as convexity and differentiability, of both and , some of which are described in the following sections. As far as the data fidelity functional is concerned, for the sake of brevity and clarity, but without restriction, we will limit the following examples to the data-fidelity functional expressed in terms of the squared .sub.2-norm as described in Equation (10).

[0057] Let us first consider gradient descent methods. In this case we consider the problem of Equation (8), where one seeks to find x∈.sup.N that minimises the cost functional (x)=(x)+λ(x). Assuming that is differentiable with respect to x, which implies that both and are differentiable with respect to x, one can express the following gradient descent iteration

x.sub.k+1←x.sub.k−γ.sub.k[∇.sub.x(x.sub.k)+λ∇.sub.x(x.sub.k)],k≥0,  (17)

where γ.sub.k∈ is a step size that may vary at each iteration k∈ ( denoting natural numbers). With suitable initialisation x.sub.0 and step size sequence (, the sequence ( may converge to the desired local or even global minimum of the cost functional F. The gradient descent algorithm consists in applying the iteration defined in Equation (17) recursively until a convergence criterion is met. The flow chart of FIG. 6 summarises the proposed image reconstruction process using the gradient descent algorithm. In step 30, the probe 3 receives or acquires a set of waveforms arising from the object of interest. In other words, in this step a set of waveforms (or signals) is received or acquired as measurements from the object of interest stimulated by (or exposed to) a physical excitation. In step 31, a measurement model is defined linking the unknown image of the medium to the recorded set of measurements. In step 32, a differentiable data-fidelity functional is constructed using H. In step 33, a differentiable regularisation functional is constructed using f.sub.θ. In step 35, the algorithm or solution is initialised as x.sub.0. In step 37, the data fidelity gradient ∇.sub.x(x.sub.k) is evaluated or determined. In step 39, the regularisation gradient ∇.sub.x(x.sub.k) is evaluated or determined. In step 41, the iteration step size γ.sub.k is evaluated or determined. In step 43, the estimate or solution is updated by computing x.sub.k+1←x.sub.k−y.sub.k [∇.sub.x(x.sub.k)+λ∇.sub.x(x.sub.k)]. In step 45, it is determined whether or not a given convergence criterion is met. In the affirmative, the process comes to an end. On the other hand, if the criterion has not been met, then the process continues in step 37.

[0058] Under strict assumptions about the cost functional , such as convexity and Lipschitz gradient, together with a suitable choice of the step size sequence (, the convergence of the gradient descent algorithm can be guaranteed. These strict assumptions concerning may be difficult to satisfy. In some situations, and/or may not be differentiable for any x∈.sup.N. In such a case, one may replace the gradient ∇.sub.x(x.sub.k) (and/or ∇.sub.x(x.sub.k)) in Equation (17) by its sub-gradient (x.sub.k)∈∂.sub.x(x.sub.k) (and/or (x.sub.k)∈∂(x.sub.k) if (and/or ) is a convex non-differentiable functional. If (and/or ) is a non-convex non-differentiable functional with local Lipschitz continuity, one may replace the gradient ∇.sub.x(x.sub.k) (and/or ∇.sub.x(x.sub.k)) in Equation (17) by its generalised gradient (x.sub.k)∈∂.sub.x(x.sub.k) (and/or (x.sub.k)∈∂.sub.x(x.sub.k)) as defined by Clarke F., (2013) Generalized gradients, “Functional Analysis, Calculus of Variations and Optimal Control”, Graduate Texts in Mathematics, vol. 264. Springer, London.

[0059] It is common to construct the data-fidelity functional as a well-behaved, convex and Lipschitz gradient functional, using a linear measurement model, but such properties may not be easily verified when constructing a regularisation functional with a trained non-linear mapping f.sub.θ, such as a NN. In this case, one may enforce the differentiability of f.sub.θ by exclusively using differentiable activation functions. If non-differentiable activation functions, such as the popular rectified linear unit (ReLU), are used and local Lipschitz continuity can be enforced, one can benefit from the generalized gradient defined above. If local Lipschitz continuity cannot be enforced, one may still approximate the gradient ∇.sub.x(x.sub.k) at discontinuous positions as it is commonly done for the training of such NNs. Having a well-guessed initialisation x.sub.0 may also be particularly advantageous. For instance, using the image obtained by the well-known DAS can provide an efficient initialisation x.sub.0 to the iteration defined in Equation (17). Many variants of the gradient descent algorithm have been developed in order to either accelerate its pace or relax the assumptions concerning for given convergence guarantees. Such variants may of course be readily applied.

[0060] The gradient descent iteration described in Equation (17) involves the gradients of both and . As an example, we consider the squared -norm data fidelity defined in Equation (10) for which the gradient can be expressed as

∇.sub.x(x.sub.k)=∇.sub.xH(x)[H(x)−y].  (18)

It can be noted that ∇.sub.xH(x) is the transpose of the Jacobian matrix defined as

[00004]$\begin{array}{cc}{J}_H\left(x\right)=\left[\frac{\partial H\left(x\right)}{\partial x_1}\frac{\partial H\left(x\right)}{\partial x_2}\mathrm{.Math.}\frac{\partial H\left(x\right)}{\partial x_N}\right]\in {ℝ}^{M\times N}.& \left(19\right)\end{array}$

[0061] In the linear case, the gradient of the squared -norm data fidelity functional becomes

∇.sub.x(x.sub.k)=H.sup.T[Hx−y].  (20)

[0062] The gradient of the regularisation functional is of particular interest. Let us consider for example the regularisation functional defined in Equation (11). Its gradient ∇.sub.x can be expressed as

∇.sub.x(x)=x−[½f.sub.θ(x)+∇.sub.xf.sub.θ(x)x],  (21)

where ∇.sub.xf.sub.θ(x)x is the directional derivative of the trained non-linear mapping f.sub.θ at position x along the direction x, where ∇.sub.xf.sub.θ(x) is the transpose of the Jacobian matrix defined as

[00005]$\begin{array}{cc}{J}_{f_{\theta }}\left(x\right)=\left[\frac{\partial f_{\theta }\left(x\right)}{\partial x_1}\frac{\partial f_{\theta }\left(x\right)}{\partial x_2}\mathrm{.Math.}\frac{\partial f_{\theta }\left(x\right)}{\partial x_N}\right]\in {ℝ}^{N\times N}.& \left(22\right)\end{array}$

[0063] However, building the Jacobian matrix J.sub.f.sub.θ(x.sub.k) at each iteration k would be computationally prohibitive. One can see that only the evaluation of the term ∇.sub.xf.sub.θ(x.sub.k)x.sub.k is required at each iteration as described in Equation (17). Provided that this evaluation is possible without having to build the Jacobian matrix, the iteration could be efficiently handled. One remarkable aspect of using a NN as f.sub.θ lies in the fact that it is differentiable in essence, or that at least means for efficiently approximating its gradient are readily available. This differentiability is a key building block of its training strategy which relies on the powerful backpropagation approach, which relies on the chain rule. Hence, one can take advantage of computationally efficient backpropagation algorithms to directly evaluate the matrix-vector product ∇.sub.xf.sub.θ(x.sub.k)x.sub.k without having to explicitly build the Jacobian matrix J.sub.fθ(x.sub.k). This strategy can be viewed as a “matrix-free” evaluation of the directional derivative and is a clear advantage when using a NN for the trained non-linear mapping f.sub.θ.

[0064] As a second example, let us consider the regularisation functional defined in Equation (14). Its gradient can be expressed as

∇.sub.x(x)=∇.sub.xr.sub.θ(x)r.sub.θ(x),  (23)

where ∇.sub.xx.sub.θ(x)r.sub.θ(x) is the directional derivative of the trained non-linear residual mapping r.sub.θ at position x along the direction r.sub.θ(x). Similarly to the previous example, one can evaluate the matrix-vector product ∇.sub.xr.sub.θ(x.sub.k)r.sub.θ(x.sub.k) at each iteration k by taking advantage of efficient backpropagation algorithms.

[0065] Projection-based methods are considered next. Let us first introduce the proximal gradient methods, which may be used to solve the problem of Equation (8) when both and are convex functionals, is differentiable and is non-differentiable. The iteration of the corresponding algorithm is given by

x.sub.k+1←pro(x.sub.k−γ.sub.k∇.sub.x(x.sub.k)),k≥0  (24)

where γ.sub.k∈ is the gradient step size at iteration k∈. Equation (24) first performs a gradient step with respect to the data fidelity ensuring closeness to the measurements and then a proximal step, which accounts for the prior knowledge about the unknown image. The proximal operator involved in Equation (24) is defined as

[00006]$\begin{array}{cc}{\mathrm{prox}}_{\mathrm{\mu ℛ}}\left(z\right)=\underset{x\in {ℝ}^N}{\mathrm{argmin}}\left\{{.Math.x-z.Math.}_2^2+\mathrm{\mu ℛ}\left(x\right)\right\},& \left(25\right)\end{array}$

for some μ∈.sup.+ and where z∈.sup.N refers to an intermediate image. Depending on the trained non-linear mapping f.sub.θ used in the construction of the regularisation functional , Problem (25) may have a closed-form solution in which case it can be readily applied. If such a closed-form solution is not available, one could solve Equation (25) iteratively which might however be computationally demanding.

[0066] Looking more closely at the proximal operator defined in Equation (25) and involved in Equation (24), it can be noted that this operation is conceptually a denoising of the image x. Taking advantage of the plug-and-play prior approach as explained for example by S. V. Venkatakrishnan, Ch. A. Bouman and Br. Wohlberg, “Plug-and-Play priors for model based reconstruction”, in 2013 IEEE Global Conference on Signal and Information Processing, 2013, pp. 945-94, one can replace the proximal operator by a suitable denoiser. Hence, replacing the proximal operator by a trained non-linear mapping f.sub.θ can provide significant improvements in the resulting image quality. The iteration defined in Equation (24) becomes

x.sub.k+1←f.sub.θ(x.sub.k−γ.sub.k∇.sub.x(x.sub.k)),k∝0.  (26)

[0067] It can be noted that in the special case where f.sub.θ is a projector onto some closed convex set , the iteration described in Equation (26) corresponds to the well-known projected gradient descent algorithm and can be used to solve the constrained optimisation problem defined in Equation (9). Hence, the iteration defined in Equation (26) may be used to form a generic projection-based algorithm, which reduces to the projected gradient descent algorithm when f.sub.θ is a projector onto some closed convex set . The generic projection-based algorithm consists in applying the iteration defined in Equation (26) recursively until a convergence criterion is met. The flow chart of FIG. 7 shows a proposed image reconstruction process using a generic projection-based algorithm. In step 50, the probe 3 receives or acquires a set of waveforms arising from the object of interest. In other words, in this step a set of waveforms (or signals) is received or acquired as measurements from the object of interest stimulated by (or exposed to) a physical excitation. In step 51, the measurement model is defined linking the unknown image of the medium to the recorded set of measurements. In step 52, is a differentiable or sub-differentiable data-fidelity functional is constructed or defined using H. In step 53, the algorithm or solution is initialised as x.sub.0. In step 55, the data fidelity gradient ∇.sub.x(x.sub.k) is evaluated or determined. In step 57, the iteration step size γ.sub.k is evaluated or determined. In step 59, the estimate or solution is updated by applying the projection x.sub.k+1←f.sub.θ(x.sub.k−γ.sub.k∇.sub.x(x.sub.k)). In step 61, it is determined whether or not a given convergence criterion is met. In the affirmative, the process comes to an end. On the other hand, if the criterion has not been met, then the process continues in step 55. When the trained non-linear mapping f.sub.θ is a general mapping it is similar to the plug-and-play approach. When the trained non-linear mapping f.sub.θ is a projector onto some closed convex set, the projection-based algorithm is equivalent to the well-known projected gradient descent algorithm.

[0068] A special case of interest arises when considering a single iteration k=0 of Equation (26) with step size γ.sub.0=1 and initialisation x.sub.0=0. Assuming a linear measurement case with a squared .sub.2-norm data fidelity functional as defined in Equation (10), for which the gradient is expressed in Equation (20), the estimate {circumflex over (x)} is given by

[00007]$\begin{array}{cc}\hat{x}=f_{\theta }\left(x_0-{\gamma }_0{H}^T\left[H{x}_0-y\right]\right)& \left(27\right)\\ =f_{\theta }\left(H^{T}y\right).& \left(28\right)\end{array}$

[0069] It has been proved that the operation H.sup.Ty is in fact the well-known DAS operation applied to the measurements y when specific assumptions are made concerning the linear measurement operator H. Hence, in this very specific case, a high-quality image estimate {circumflex over (x)} can be obtained from a low-quality image estimate {tilde over (x)} reconstructed by DAS such that {tilde over (x)}=H.sup.T y by a simple direct mapping

{circumflex over (x)}=f.sub.θ({tilde over (x)}).  (29)

[0070] Thus, the low-quality image estimate {tilde over (x)}=H.sup.Ty is first computed and then the estimate is computed by applying the mapping, such that {circumflex over (x)}=f.sub.θ({tilde over (x)}).

[0071] Such a simplification has the great advantage of only requiring the reconstruction of a low-quality image by means of DAS followed by a single evaluation of f.sub.θ. Hence, the trained non-linear mapping f.sub.θ is approximating the inverse of the mapping F defined in Equation (3). Even though simplistic, this special case may already provide significant improvements on the reconstructed high-quality image {circumflex over (x)}, which is significantly enhanced compared to the low-quality image {tilde over (x)}.

[0072] It can be noted that both image reconstruction examples presented above (gradient descent and projection-based methods) rely at each iteration of an iterative algorithm on a globally trained non-linear mapping f.sub.θ. However, one can readily specialise this approach with a set of non-linear mappings

[00008]${\left\{f_{{\theta }_k}\right\}}_{k\in \mathbb{N}}$

in which each f.sub.θ.sub.k is specifically designed and trained for the given level of artefacts present at each iteration k∈ of an iterative algorithm. In this case the image reconstruction method would be based on at least two trained non-linear mappings (a multiplicity of trained non-linear mappings).

[0073] Above it was described a general method for solving an optimisation problem, such as the one described in Equation (8), where the regularisation functional is constructed using a trained non-linear mapping f.sub.θ. In the following, some insights are provided on how to design and train such a non-linear mapping in order to fully benefit from the proposed approach.

[0074] As a non-limiting example, the residual “encoder-decoder” architecture shown in FIG. 5 is considered. Such an architecture may be implemented as a deep residual convolutional NN, such as the one presented in FIG. 8. The encoding of an input image x is performed through the downward path of the NN architecture (from top-left to bottom-centre) and the decoding is performed in the upward path of the NN architecture (from bottom-centre to top-right) up to the intermediate residual output r.sub.θ(x). The residual output is then summed to the input image to form the final output f.sub.θ(x)=x+r.sub.θ(x). The solid arrows specify operations, while the dashed arrows specify intrinsic skip connections. The dashed-dotted arrows do not represent an operation but represent a transfer of quantity. Each of the dotted boxes frames a decomposition level. In this specific example, three decomposition levels are considered, namely l=0, 1, 2. The rectangles with solid outlines represent feature maps, i.e. the output of an operation or layer, at the respective depths of the NN. The width (i.e. the dimension along a horizontal line) of each box indicates the channel number present in the feature map. The first operation expands the channel number of the input image x to a user-definable channel number N.sub.c using a convolutional operation with a kernel size of (1×1). The last operation before the residual output r.sub.θ(x) performs a contraction of the channel number N.sub.c back to the channel number of the input image using a convolutional operation with kernel size of (1×1). Within each decomposition level, a (3×3) convolutional operation followed by a non-linear activation (nl. in FIG. 8) is performed between each feature map without changing the dimension of the feature map. Such a non-linear activation may for example be the well-known rectified linear unit (ReLU), the exponential linear unit (ELU) or the sigmoid function. Between each decomposition level in the downward path, a downsampling operation is performed using a multiscale (3×3) convolutional operation with (or without) an expansion of the number of channels, followed by a non-linear activation. Such a multiscale operation may for example be the combination of a max pooling operation followed by a convolution operation. It may also be a strided convolution with (or without) an expansion of the number of channels. Between each decomposition level in the upward path, an upsampling operation is performed using a multiscale (3×3) convolutional operation followed by a non-linear activation. Such a multiscale operation may for example be a transposed and strided convolution operation with (or without) a contraction of the number of channels.

[0075] All the filter coefficients of the NN architecture are trainable parameters. These parameters, also called weights, may be trained by minimising the empirical risk

[00009]$\begin{array}{cc}\hat{\theta }=\underset{\theta }{\mathrm{argmin}}\frac{1}{L}{.Math.}_{j=1}^Lℒ\left(f_{\theta }\left({\tilde{x}}_j\right),{\hat{x}}_j\right),& \left(30\right)\end{array}$

where is a training loss function, {{circumflex over (x)}.sub.j}.sub.j=1.sup.L and {{tilde over (x)}.sub.j}.sub.j=1.sup.L are sets of L images. It is common to construct the sets {{circumflex over (x)}.sub.j}.sub.j=1.sup.L and {{tilde over (x)}.sub.j}.sub.j=1.sup.L by respectively sampling from a high-quality subspace V and a low-quality subspace W. Such samples may for example be acquired experimentally or by numerical simulation. The set {{circumflex over (x)}.sub.j,{tilde over (x)}.sub.j}.sub.j=1.sup.L is often referred to as a training set. Optimisation algorithms, such as stochastic gradient descend (SGD), batch SGD or Adam optimizer (D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization”, pp. 1-15, 2014. arXiv: 1412.6980), which is an SGD variant using momentum evaluated from previous iterations, may for instance be used for solving the minimisation problem defined in Equation (30). To this end, the gradient of the loss function with respect to the trainable parameters ∇.sub.θ may for example be computed by backpropagation using the popular chain rule. A typical loss function is the mean squared error (MSE) defined as

(f.sub.θ({tilde over (x)}.sub.j),{circumflex over (x)}.sub.j)=½∥{circumflex over (x)}.sub.j−f.sub.θ({tilde over (x)}.sub.j)∥.sub.2.sup.2.  (31)

[0076] However, different loss function, risk minimisation function and optimisation algorithm may be used.

[0077] For the non-linear mapping f.sub.θ, a neural network architecture called U-Net may be used (among other possibilities) which was first introduced by O. Ronneberger, Ph. Fischer, and Th. Brox, “U-Net: Convolutional Networks for Biomedical Image Segmentation”, pp. 1-8, 2015. arXiv: 1505.04597 for biomedical image segmentation purpose and recently adapted in the context of sparse-view CT image reconstruction. It is in essence very similar to the residual non-linear mapping presented in FIG. 8. The parameters θ of f.sub.θ are optimised using the Adam optimiser while the MSE, which was described earlier, is used as a training loss.

[0078] The training set may be composed of approximatively 1000 to 100 000 or 10 000 to 30 000 images generated using a specific simulator for example. Synthetic phantoms may be used for this purpose and may be randomly generated with inclusions of variable size and echogenicity. The probability distributions for inclusion size and echogenicity may be modelled to closely mimic the distribution present in the object of interest. Numerical phantoms may be designed and used to obtain a fully diffusive speckle pattern throughout the entire image with an echogenicity ranging from −∞ dB (anechoic regions) to 30 dB (hyperechoic regions). The main advantages of using simulated data rather than experimental data are as follows. It provides full control over the distribution of artefacts and inclusions of different echogenicity. This is useful for training a neural network because, if a specific artefact is only scarcely represented in the training set, the neural network will put very little focus in removing it. It is also possible to impose a wider diversity in the training set than what would be possible with experimental data. The fact that the simulation environment is entirely controlled also prevents blurring effects induced by tissue motion. Another advantage is the possibility to simulate non-existing US imaging configurations. For example, the LQ configuration defined in Table 1 is based on a realistic probe, but the probe used in the HQ configuration does not exist and allows the reconstruction of ideal images for the given aperture. The main disadvantage of using simulated training set is that it is not self-evident whether or not the results will translate when applied to experimental data.

[0079] Once trained, the direct mapping approach given in Equation (29) may for instance be used. A high-quality image estimate {circumflex over (x)} can thus be obtained by simply applying the trained non-linear mapping f.sub.θ to a low-quality image estimate {tilde over (x)}.

[0080] To summarise, the present invention is motivated by the fact that reconstructing images from a high-quality US imaging configuration, in an ultrafast US imaging context, comes at a price and might not even be feasible in practice. Indeed, high-quality ultrafast images are generally obtained by averaging or compounding multiple images of lower quality reconstructed sequentially. An example of this strategy is the well-known coherent PW compounding where steered PWs are transmitted sequentially and the final image is obtained by averaging all images reconstructed for each steered PW. Since transmit-receive events are sequential, the round-trip time-of-flight of the acoustic wave can already be a limiting factor in scenarios where very high frame rates are required. The computational requirements for reconstructing many images to be compounded may also not be viable in some situations such as in 3D imaging or low-power environments. Moreover, the motion of the object (e.g. a human tissue) can induce a substantial blurring of the resulting image when averaging sequential images reconstructed at different time instants. Thus, one of the goals of the present invention is to reconstruct high-quality images with a single or very few PW measurements.

[0081] While the invention has been illustrated and described in detail in the drawings and foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive, the invention being not limited to the disclosed embodiment. Other embodiments and variants are understood, and can be achieved by those skilled in the art when carrying out the claimed invention, based on a study of the drawings, the disclosure and the appended claims. For example, the teachings of the present invention can be applied to photoacoustic imaging, where the transmitted acoustic waves are replaced with electromagnetic-radiation pulses.

[0082] In the claims, the word “comprising” does not exclude other elements or steps, and the indefinite article “a” or “an” does not exclude a plurality. The mere fact that different features are recited in mutually different dependent claims does not indicate that a combination of these features cannot be advantageously used.