Model-based image reconstruction method

Abstract

The present invention concerns a pulse wave image reconstruction method to be used for example in ultrasound imaging. The proposed method is based on an image measurement model and its adjoint operator. The proposed method introduces matrix-free formulations of the measurement model and its adjoint operator. The proposed method has the advantage that the reconstructed image has a very high quality and that it can be reconstructed quickly.

Description

TECHNICAL FIELD

The present invention relates to a model-based image reconstruction method suitable to be used for example in ultrasound imaging. The invention also relates to an imaging apparatus for carrying out the method.

BACKGROUND OF THE INVENTION

Pulse wave imaging applications, such as ultrasound (US) imaging applications, may be used in various technical fields. In medical applications, it is used to see internal body structures such as tendons, muscles, joints, vessels and internal organs. However, pulse wave imaging techniques are also widely used outside the medical field. One commercial use of ultrasound is non-destructive testing (NDT), which is a process of inspecting, testing, or evaluating materials, components or assemblies for discontinuities, or differences in characteristics without destroying the serviceability of the part or system. In other words, when the inspection or test is completed, the part can still be used. An example of an NDT would be in oil pipes. In this case, ultrasound waves are used to propagate down the pipe, and in this manner, it is possible to detect if there is a crack or defect as a reflection of the pulse. With knowledge of the speed of the wave in the pipe and the time interval between the reflection from the crack and sending the waves, the position of the defect can be estimated. The same theory used for crack detection in pipes can be used to detect cracks in aircraft structures/wing panels etc and to detect the build-up of ice on the surface of the wing for example. There are of course many other possible NDT applications for the pulse wave imaging. Pulse wave imaging has gained popularity, because of its safety, portability and real-time capability.

Classical US imaging uses multiple focused time-separated transmit beams to reconstruct images along scan lines. The frame rate is therefore limited by the number of focused transmit beams and cannot exceed a few tens of images per second. Although this frame rate is sufficient for many applications, the understanding of more complex phenomena, such as the propagation of shear waves for elastography, requires higher frame rates.

One way to overcome this limitation consists in decreasing as much as possible the number of transmit beams. Synthetic aperture (SA) methods try to overcome this problem by using only few transducer elements to sequentially insonify the whole medium. Another option is ultrafast US (UFUS) imaging, which exploits the idea of using plane waves (PW) or diverging waves (DW) to insonify the whole field-of-view at once, allowing US systems to reach thousands of frames per second.

An important limitation of such approaches is a degraded image quality. Indeed, compared to focused transmit beams where the energy is concentrated in a specific region of interest, the energy of a PW or a DW is spread over the insonified medium, resulting in lower signal-to-noise ratio and poor spatial resolution. One way to address this problem consists in averaging multiple images, obtained with different insonification angles in the context of PW imaging and with different point source locations in the context of DW imaging for UFUS or with different groups of transducer elements for SA imaging. In UFUS, this process is called coherent compounding. While the implementation of such a technique is rather simple, it requires multiple insonifications thus reducing the frame rate.

An alternative to averaging multiple insonifications consists in using more efficient reconstruction methods than classical techniques, which revolve mainly around delay-and-sum (DAS) beamforming. One popular group of methods relies on the use of iterative algorithms to solve the ill-posed image formation problem induced by US image reconstruction. These methods are built upon forward models of the problem. The main problem with the known solutions resides in their computational complexity, usually translated into storage requirements of the corresponding matrix representation. Some models proposed require the storage of several hundreds of giga bytes for matrix coefficients in two dimensions (2D). This issue severely limits the appeal of the iterative methods in view of the classical approaches.

SUMMARY OF THE INVENTION

It is an object of the present invention to overcome at least some of the problems identified above related to image reconstruction in pulse wave imaging.

According to a first aspect of the invention, there is provided an image reconstruction method as recited in claim 1.

The proposed method has the advantage that the reconstructed image has a very high quality and that it can be reconstructed quickly. In other words, compared to prior art solutions, image contrast and resolution can be significantly improved.

According to a second aspect of the invention, there is provided an imaging apparatus configured to carry out the method according to the first aspect of the invention as recited in claim 15.

Other aspects of the invention are recited in the dependent claims attached hereto.

BRIEF DESCRIPTION OF THE DRAWINGS

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

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;

FIG. 2 is a schematic illustration showing the imaging configuration used in the proposed imaging method according to an example of the present invention;

FIG. 3 is a schematic illustration showing the imaging configuration and a continuous domain reflectivity hypersurface used in the proposed imaging method according to an example of the present invention;

FIG. 4 is a schematic illustration showing the imaging configuration and a discrete domain reflectivity hypersurface used in the proposed imaging method according to an example of the present invention;

FIG. 5 is a flow chart summarising the image reconstruction method according to an example of the present invention;

FIG. 6 is a flow chart illustrating the image reconstruction method in more detail according to an example of the present invention;

FIG. 7 is a flow chart illustrating the image reconstruction method in more detail according to another example of the present invention;

FIG. 8 is a flow chart illustrating a process of implementing an image measurement model according to an example of the present invention; and

FIG. 9 is a flow chart illustrating a process of implementing an adjoint operator of the measurement model according to an example of the present invention.

DETAILED DESCRIPTION OF AN EMBODIMENT OF THE INVENTION

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 US imaging, but the teachings of the invention are not limited to this environment. The teachings of the invention are equally applicable in other fields where pulse wave imaging can be used. Identical or corresponding functional and structural elements which appear in the different drawings are assigned the same reference numerals.

The present invention in the present embodiment proposes an image formation method and apparatus or system in the context of 2D and 3D ultrasound imaging based on solving an ill-posed linear inverse problem of the following form:m=Hγ+n, where m are the measurements, γ is the object under scrutiny or more specifically its reflectivity function, H is the measurement model and n is the noise. The proposed image formation method relies on two main pillars: a fast and matrix-free measurement model H and an image reconstruction method which permits to retrieve an estimate of the object under scrutiny γ given the measurements m. The object may be an internal part or element of a larger structure. The proposed method can be carried out by a standard US system comprising a US probe, an image formation module and a post-processing and display modules, as described in FIG. 1.

FIG. 1 is 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 system 1 comprises an ultrasound probe 3, which in this example comprises a linear array of piezoelectric elements or sensors, 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.

The US probe 3 is connected to an image formation module or apparatus 5, which is configured to carry out the image formation. A measurement model estimated in a measurement model unit 7 is used in the image formation process. The proposed measurement model as explained later in detail is based on a pulse-echo spatial impulse response model. The present invention introduces matrix-free formulations of the measurement model and its adjoint, the adjoint of an operator being defined as a continuous extension of the transpose of a square matrix in the case of a continuous operator, which represent a novel image formation method in pulse wave imaging. The description of the measurement model may be achieved through two main steps:

• • Parametric equations of the round-trip time-of-flight of the pulse waves are derived. It will be shown that the model can be recast as an integral over hypersurfaces that are specified as explained later.
  • The proposed integral may then be discretised as described later.

The model may be derived for example for a PW, DW and SA imaging. It is compatible with radio frequency and in-phase and quadrature (IQ) data. Finally, a matrix-free formulation of the adjoint of the proposed measurement model is derived. It will be demonstrated that it can also be recast as an integral over hypersurfaces, defined as 1D manifolds in the case of 2D imaging and as 2D manifolds in the case of 3D imaging, that are parametrised and that the same discretisation scheme as proposed for the measurement model may be used.

The image formation apparatus 5 also comprises an image reconstruction unit 9 for running an image reconstruction process or method. The image reconstruction method aims to retrieve an estimate {circumflex over (γ)} of the image under scrutiny. The estimate can be written as:{circumflex over (γ)}=f(m,H,Ξ),where f accounts for the image reconstruction process with hyperparameters Ξ. Two different image reconstruction methods are presented in the present description:

• • Analytical approaches where f(m, H, Ξ) has an explicit formulation (eg backpropagation methods using the adjoint operator of H).

Iterative approaches where f(m, H, Ξ) is an iterative algorithm, which solves the following problem:

$\hat{\gamma }=\underset{\gamma }{\mathrm{argmin}}𝒥\left(\gamma ,m\right)=ℱ\left(H\gamma ,m\right)+\mathrm{\lambda ℛ}\left(\gamma \right),$

• • where (γ, m) denotes an objective function which comprises a functional , also referred to as a data discrepancy term, measuring the error between Hγ and m, denotes a non-negative functional (regularisation) and λ denotes a regularisation parameter. In this case, Ξ are the hyperparameters related to the optimisation algorithm used to solve the minimisation problem, eg the number of iterations, the stopping criterion of the algorithm etc.

After image reconstruction, the reconstructed data are sent to a post-processing and display unit or module 11. The post-processing step covers a range of applications, such as B-mode imaging, colour Doppler imaging, vector Doppler imaging and elastography.

For B-mode imaging using RF data, an envelope detection is applied onto the reconstructed data. 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. If B-mode imaging is used with IQ data, the magnitude of the signal is extracted. The envelope detection step is followed by a normalisation and dynamics compression steps. For Doppler imaging and elastography, the reconstructed RF or IQ data are directly used without post-processing.

Before explaining the proposed method in detail, some notations used in the present description are explained first. In the pulse-echo configuration considered in the present description, echoes from a medium Ω∈³, characterised by its reflectivity (function) γ(r) with r=[x, y, z]^T∈Ω are detected by sensors 13. FIG. 2 illustrates a specific simplified measurement configuration for PW imaging. The following notations are used throughout the present description:

Mathematical or Similar Notations:

• • Vectors are denoted with a bold font.

${v}_p=\sqrt[p]{{\Sigma }_{i=1}^N|v_i{|}^p}$ denotes the -norm of a vector v∈^N, where denotes the space of real numbers.

• • The space of the square integrable functions of a Hilbert space W is denoted as L₂(W). The inner product between two functions ƒ, g∈L₂(W) is defined as ƒ|g_w=∫_x∈Wƒ(x)g(x)dx.

Ultrasound Related Notations:

• • m(p, t) is the electrical signal received at a time instant t by the sensor 13 positioned at p. In other words, m(p, t) defines one sample of a continuous set of measurements (equivalently one sample of element raw data), computed at a measurement grid point (p, t).
  • γ(r) defines one value of a set of reflectivity values (equivalently image samples) computed at a reflectivity grid point r.
  • Ω denotes the medium or object of interest, ie the continuous set of locations of the reflectivity samples.
  • o(r, p) is the element directivity between the sensor located at p and the point located at r.
  • a(r, p) defines the value of the attenuation between a point source located at p and a point located at r.
  • v_pe(t) is the received pulse shape (ie the pulse-echo waveform).
  • t_Tx(r, q) denotes the transmission propagation delay, ie the time taken by the pulse wave transmitted by a transmitting device or element 13 located at q to reach the point located at r.
  • t_Rx(r, p) denotes the receive propagation delay, ie the time taken by the pulse-echo wave transmitted by a point located at r to reach the sensor located at p.

Measurement Grid Related Notations:

• • A set of transmitting elements: (X^j)_j=1^Nelt, where N_elt refers to the number of transmitting elements, also referred to as transmitters.
  • A set of receive sensors: (Πⁱ)_i=1^Nelr, where N_elr refers to the number of receive sensors, also referred to as receivers.
  • A set of point source locations inside a transmit element: X_d^j={q^k∈X^j, k={1 . . . N_q^j}}.
  • A set of point source locations inside a receive sensor: Π_dⁱ={pⁿ∈Πⁱ, n={1 . . . N_pⁱ}}.
  • A set of time samples: T_d={t^l∈, l={1 . . . N_t}}.

Reflectivity Grid Related Notations:

• • A set of locations of reflectivity values: Ω_d={rⁿ∈Ω, n={1 . . . N_γ}} where N_γ refers to the number of reflectivity values.

The proposed method defines a first grid and a second grid which both lie in the 3D space. The first grid is referred to as a reflectivity grid, which comprises reflectivity grid points at intersections of two crossing grid lines, which in this example form an angle of substantially 90 degrees with respect to each other. Every reflectivity grid point of the reflectivity grid belongs to the medium Ω. The second grid is referred to as a measurement grid, which comprises measurement grid points at intersections of two crossing grid lines, which in this example form an angle of substantially 90 degrees with respect to each other. The measurement grid points coincide with the locations of the receiving sensors, in one dimension, and with the sampling time instants of each sensor in the other dimension.

Estimation, evaluation, computation or implementation of the measurement model at a given measurement grid point according to the present invention is explained next. In this example, the measurement model is estimated at each measurement grid point. The measurement model used in this disclosure is a pulse-echo spatial impulse response model first introduced by G. E. Tupholme, “Generation of acoustic pulses by baffled plane pistons”, Mathematika, vol. 16, p. 209, 1969, and P. R. Stepanishen, “The Time-Dependent Force and Radiation Impedance on a Piston in a Rigid Infinite Planar Baffle”, J. Acoust. Soc. Am., vol. 49, p. 76, 1971. According to this model, one can write the below equations which account for the raw data sample recorded at position p and time t.

A continuous parametric formulation of the measurement model for a configuration of one transmitting point source element and one receiving point source sensor is explained next.

The measurement model is defined as:

$m\left(p,t\right)={\int }_{r\in \Omega }\frac{o\left(r,p\right)o\left(r,q\right)a\left(r,p\right)a\left(r,q\right)}{4{\pi }^2{r-p}_2{r-q}_2}\gamma \left(r\right)v_{pe}\left(t-t_{Tx}\left(r,q\right)-t_{Rx}\left(r,p\right)\right)\mathrm{dr}.$ The above equation can be expressed as shown below, and defines the measurement model as an integral along the hypersurface Γ(p, q, t):

$m\left(p,t\right)=\int {\int }_{\tau \in ℝ,r\in \Gamma \left(p,q,\tau \right)}\frac{o_d\left(r,p,q\right)\gamma \left(r\right)}{\left|{\nabla }_rℊ\right|}{v}_{pe}\left(t-\tau \right)d\sigma \left(r\right)d\tau ,\mathrm{where}$ $o_d\left(r,p,q\right)=\frac{o\left(r,p\right)o\left(r,q\right)a\left(r,p\right)a\left(r,q\right)}{4{\pi }^2{r-p}_2{r-q}_2}$ $\Gamma \left(p,q,t\right)=\left\{r\in \Omega |ℊ\left(r,p,q,t\right)=0\right\}$ $ℊ\left(r,p,q,t\right)=t-t_{Tx}\left(r,q\right)-t_{Rx}\left(r,p\right)\left|{\nabla }_rℊ\right|=\left|\frac{\partial ℊ}{\partial r}\right|$

In the above equations, the hypersurface, or more specifically a reflectivity hypersurface, is Γ(p, q, t), defined by the equation g(r, p, q, t)=0, where the variable of interest is r. Γ(p, q, t) defines a conic when r lies in the 2D plane (parabola or ellipse) and a quadric when r lies in a 3D volume. The same shapes apply also to the measurement hypersurface. FIG. 3 schematically illustrates the measurement configuration in continuous domain showing also the reflectivity hypersurface.

The reflectivity hypersurface Γ(p, q, t) can now be parametrised in order to reduce the integration dimension. The reflectivity hypersurface Γ(p, q, t), defined by the equation g(r, p, q, t)=0, can equivalently be defined by a set of reflectivity parametric equations. In mathematics, parametric equations define a set of quantities as functions of one or more independent variables (ie parameters). Parametric equations are commonly used to express the coordinates of the points that make up a geometric object such as a surface or curve, in which case the equations are collectively called a parametric representation or parameterisation of the object. The set of reflectivity parametric equations can now be defined as follows:r∈Γ(p,q,t)⇔r=r(α,p,q,t).

The equation on the right of the equivalence defines the reflectivity parametric equations of r since r now depends on a parameter α∈A, where A⊂ if r∈³ and A⊂ when r∈². It is to be noted if there is one set of reflectivity parametric equations per measurement value or sample. Furthermore, the hypersurface can be obtained by drawing the set of reflectivity parametric equations. Equipped with the above formulation of the set of reflectivity parametric equations, it is possible to derive the continuous parametric formulation of the measurement model H as follows:

$m\left(p,t\right)=\int {\int }_{\tau \in ℝ,\alpha \in A}{o}_d\left(r\left(\alpha ,p,q,\tau \right),p,q\right)\gamma \left(r\left(\alpha ,p,q,\tau \right)\right)J\left(\alpha ,p,q,\tau \right)v_{\mathrm{pe}}\left(t-\tau \right)d\alpha d\tau ,\mathrm{where}J\left(\alpha ,p,q,\tau \right)=\frac{{\nabla }_{\alpha }r}{{\nabla }_rℊ}{\nabla }_{\alpha }r=\frac{\partial r}{\partial \alpha }.$

A continuous measurement model for a configuration of one transmitting finite length element of arbitrary shape and one receiving point source sensor is explained next. In this case, the transmitting element has a finite length and an arbitrary shape. It is characterised by a continuous set of point source locations X^j. The measurement model is defined asm(p,t)=∫∫_{q∈Xj,r∈Ω}o_d(r,p,q)γ(r)v_pe(t−t_Tx(r,q)−t_Rx(r,p))drdσ(q).

It can be seen that the main difference with the model described in connection with the configuration of one transmitting point source element and one receiving point source sensor lies in the fact that there is one more integral, which corresponds to the sum of the contributions of all the point sources which form the finite length sensor of arbitrary shape.

The corresponding parametric formulation is given by:m(p,t)=o_d(r(α,p,q,τ),p,q)γ(r(α,p,q,τ))|J(α,p,q,τ)|v_pe(t−τ)dαdσ(q)dτ.

A continuous measurement model for a configuration of one transmitting finite length element of arbitrary shape and one receiving finite length sensor is explained next. In this case, the receiving sensor has a finite length and arbitrary shape. It is defined by the set of point source locations Πⁱ. The coordinate of the resulting signal is denoted by m(ξⁱ, t), where ξⁱ∈Πⁱ. It can be for example the coordinate of the point in the centre of Πⁱ.

It can be written: m(ξⁱ, t)=∫_p∈Πim(p, t)dσ(p). The resulting measurements are obtained by integrating along the transducer surface defined by Πⁱ. When integrating the equation which defines m(p,t) into the above integral, the following relationship holds:m(ξⁱ,t)=∫∫∫_{p∈Πi,q∈Xj,r∈Ω}o_d(r,p,q)γ(r)v_pe(t−t_Tx(r,q)−t_Rx(r,p))drdσ(q)dσ(p).

The corresponding parametric formulation is given by:m(ξⁱ,t)=o_d(r(α,p,q,τ),p,q)γ(r(α,p,q,τ))|J(α,p,q,τ)|v_pe(t−τ)dαdσ(q)dσ(p)dτ.

A continuous measurement model for a configuration of multiple transmitting finite length elements of arbitrary shape and one receiving finite length sensor of arbitrary shape is explained next. In this case, a set of point source locations X^j is associated with the transmitting element j, leading to the following equation:m(ξⁱ,t)=Σ_j=1^Nelt∫∫∫_{p∈Πi,q∈Xj,r∈Ω}o_d(r,p,q)γ(r)v_pe(t−t_Tx(r,q)−t_Rx(r,p))drdσ(q)dσ(p).

This equation is highly similar to the one described in connection with the configuration of one transmitting finite length elements of arbitrary shape and one receiving finite length sensor of arbitrary shape. The only difference is in the sum that has been introduced to take into account the multiple transmitting elements.

The corresponding parametric formulation is given by:m(ξⁱ,t)=Σ_j=1^Nelto_d(r(α,p,q,τ),p,q)γ(r(α,p,q,τ))|J(α,p,q,τ)|v_pe(t−τ)dαdσ(q)dσ(p)dτ.

A continuous measurement model for a configuration of multiple transmitting finite length elements of arbitrary shape and multiple receiving finite length sensors of arbitrary shape is explained next.

When a set of multiple receiving sensors is used, a set of point source locations Πⁱ and a coordinate ξⁱ is associated with the receiving sensor i. Thus, the receiving measurement field can be defined asm(ξ,t)=Σ_i=1^Nelrm(ξ,t)δ(ξ−ξⁱ),which is only defined for discrete coordinates ξⁱ. In the remainder of the present description, the following holds: ξ∈P⊂³.

The continuous parametric formulation of the measurement model is next discretised. This step consists in discretising the parametric formulation of the measurement model described in connection with the configuration of multiple transmitting finite length elements of arbitrary shape and one receiving finite length sensor of arbitrary shape. However, any of the models described above could be discretised instead. By discretisation, it is meant that we select some discrete values α^p of the parameter α, which defines the following set of parameters: A_d={α^p∈A, p={1 . . . N_α}}. It can be seen that the selection of the discrete values of the parameter α leads to a discretisation of the hypersurface Γ(p, q, t), which will now be described by only N_α points (r(α^p, p, q, t))_αp∈Ad. FIG. 4 schematically illustrates the measurement configuration in discrete domain showing also the reflectivity hypersurface. The discretisation of the parametric formulation of the measurement model can be expressed as:m(ξⁱ,t^l)=({tilde over (m)}(ξⁱ)*v_pe)(t^l),where v_pe∈^Nt is the pulse shape evaluated at each time sample of the measurement grid and {tilde over (m)}(ξⁱ)=({tilde over (m)}(ξⁱ, t^l))_tl∈Td with

$\tilde{m}\left({\xi }^i,t^l\right)=\sum _{j=1}^{N_{\mathrm{elt}}}\sum _{q^k\in X_d^j}\sum _{p^n\in {\Pi }_d^i}\sum _{{\alpha }^p\in A_d^j}{o}_d\left(r\left({\alpha }^p,p^n,q^k,t^l\right),p^n,q^k\right)\gamma \left(r\left({\alpha }^p,p^n,q^k,t^l\right)\right)J\left({\alpha }^p,p^n,q^k,t^l\right)w\left(q^k\right)u\left({\alpha }^p\right),$ where w(q^k), u(α^p) and z(pⁿ) are integration weights (related to the discretisation of the continuous integrals). It can be noticed that, for each point q^k of the transmitting element, the above equation comprises a set of reflectivity values (γ(r(α^p, pⁿ, q^kt^l)))_αp∈Adj, referred to as hypersurface reflectivity samples or values, which are estimated on the points of the reflectivity hypersurface.

The hypersurface reflectivity values are next interpolated to substantially coincide with the reflectivity grid points. The hypersurface reflectivity values are usually unknown since they do not lie on the reflectivity grid points. Thus, they need to be approximated from the reflectivity samples defined on the reflectivity grid by means of an interpolation kernel as introduced by P. Thévenaz, T. Blu and M. Unser, “Interpolation revisited”, IEEE Transactions on Medical Imaging, vol. 19, p. 739-758, 2000. Let us thus introduce the interpolation kernel φ:Ω→ such that:γ(r(α^p,pⁿ,q^k,t^l))≈Σ_s=1^Nγγ^sφ(r^s−r(α^p,pⁿ,q^k,t^l)).

The above interpolation equation relates the hypersurface reflectivity values γ(r(α^p, pⁿ, q^k, t^l)) (which are typically not on the reflectivity grid points) to the reflectivity samples γ^s.

A summation of the discretised and interpolated hypersurface reflectivity samples or values is carried out next. The interpolation equation is inserted into the discretised measurement model and the following summation is performed:{tilde over (m)}(ξⁱt^l)=Σ_j=1^NeltΣ_qk∈XdjΣ_pn∈ΠdiΣ_αp∈Adjo_d(r(α^p,pⁿ,q^k,t^l),pⁿ,q^k)(Σ_s=1^Nγγ^sφ(r^s−r(α^p,pⁿ,q^k,t^l)))|J(α^p,pⁿ,q^k,t^l)|w(q^k)z(pⁿ)u(α^p).

In this example, the output of the summation is eventually convolved with the pulse shape, as described in the previous equations in order to obtain an estimate of the measurement value, also referred to as a measurement value estimate, at the considered measurement grid point. The above described process is carried out for every measurement value. Thus, the outcome of the estimation or processing of the measurement model is the estimate of the measurement value at the respective measurement grid point.

The estimation, evaluation, computation or implementation of the adjoint operator of the measurement model at each reflectivity grid point is carried out next. First, a continuous parametric formulation of the adjoint operator of the measurement model is obtained. The continuous formulation of the adjoint operator can be derived from the measurement model, ie from example from any of the continuous measurement models given above, using functional analysis tools. The equations below are derived from the case of multiple finite length transmitting elements of arbitrary shape and multiple finite length receiving sensors of arbitrary shape. Indeed, if we define an operator : L₂(Ω)→L₂(P×T) such that {γ}(ξ, t)=Σ_i=1^Nelrm(ξ, t)δ(ξ−ξⁱ), then the adjoint operator : L₂(P×T)→L₂(Ω) of is defined as:n|{γ}_P×T={n}|γ_Ω,∀n∈L₂(P×T),∀γ∈L₂(Ω).

Let us consider the following inner product:

${\langle n|\left\{\gamma \right\}\rangle }_{P\times T}=\int {\int }_{\xi \in P,t\in T}n\left(\xi ,t\right)\left\{\gamma \right\}\left(\xi ,t\right)d\xi \mathrm{dt}$ ${\langle n|\left\{\gamma \right\}\rangle }_{P\times T}=\int {\int }_{\xi \in P,t\in T}n\left(\xi ,t\right)\sum _{i=1}^{N_{\mathrm{elr}}}m\left(\xi ,t\right)\delta \left(\xi -{\xi }^i\right)d\xi \mathrm{dt}$ ${\langle n|\left\{\gamma \right\}\rangle }_{P\times T}=\sum _{i=1}^{N_{\mathrm{el}}}{\int }_{t\in T}n\left({\xi }^i,t\right)m\left({\xi }^i,t\right)\mathrm{dt}$ ${\langle n|\left\{\gamma \right\}\rangle }_{P\times T}=\sum _{i=1}^{N_{\mathrm{elr}}}{\int }_{t\in T}n\left({\xi }^i,t\right)\left[\sum _{j=1}^{N_{\mathrm{elt}}}\int \int {\int }_{q\in X^j,p\in {\Pi }^i,r\in \Omega }^{}{o}_d\left(r,p,q\right)\gamma \left(r\right)v_{\mathrm{pe}}(t-t_{\mathrm{Tx}}\left(r,q\right)-t_{\mathrm{Rx}}\left(r,p\right))\mathrm{drd}\sigma \left(q\right)d\sigma \left(p\right)\right]\mathrm{dt}$ ${\langle n|\left\{\gamma \right\}\rangle }_{P\times T}={\int }_{r\in \Omega }\gamma \left(r\right)\sum _{j=1}^{N_{\mathrm{elt}}}\sum _{i=1}^{N_{\mathrm{elr}}}\int {\int }_{q\in X^j,p\in {\Pi }^i}{o}_d\left(r,p,q\right)[{\int }_{t\in T}^{}n\left({\xi }^i,t\right)v_{\mathrm{pe}}\left(t-t_{\mathrm{Tx}}\left(r,q\right)-t_{\mathrm{Rx}}\left(r,p\right)\right)\mathrm{dt}]d\sigma \left(q\right)d\sigma \left(p\right)\mathrm{dr}$ ${\langle n|\left\{\gamma \right\}\rangle }_{P\times T}={\int }_{r\in \Omega }\gamma \left(r\right)\left[\sum _{j=1}^{N_{\mathrm{elt}}}\sum _{i=1}^{N_{\mathrm{elr}}}\int {\int }_{q\in X^j,p\in {\Pi }^i}{o}_d\left(r,p,q\right)\left(n\left({\xi }^i\right)*u\right)(t_{\mathrm{Tx}}\left(r,q\right)+t_{\mathrm{Rx}}\left(r,p\right))d\sigma \left(q\right)d\sigma \left(p\right)\right]\mathrm{dr},$ where u(t)=v_pe(−t) is the matched filter of the pulse echo waveform. We now have:n|{γ}_P×T=∫_r∈Ωγ(r){n}(r)dr={n}|γ_Ω,where{n}(r)=Σ_j=1^NeltΣ_i=1^Nelr∫∫_{q∈Xj,p∈Πi}o_d(r,p,q)(n(ξⁱ)*u)(t_Tx(r,q)+t_Rx(r,p))dσ(q)dσ(p).

Equipped with the adjoint operator of the measurement model, we define the estimate of the reflectivity value at the reflectivity grid point r, given an estimate of the measurement value, as:

$\hat{\gamma }\left(r\right){†}^{}\left\{m\right\}\left(r\right)$ $\hat{\gamma }\left(r\right)=\sum _{j=1}^{N_{\mathrm{elt}}}\sum _{i=1}^{N_{\mathrm{elr}}}\int {\int }_{q\in X^j,p\in {\Pi }^i}{o}_d\left(r,p,q\right)\left(m\left({\xi }^i\right)*u\right)\left(t_{\mathrm{Tx}}\left(r,q\right)+t_{\mathrm{Rx}}\left(r,p\right)\right)d\sigma \left(q\right)d\sigma \left(p\right).$

This case is simpler than the measurement model. Indeed, the parametric formulation is direct and the measurement parametric equation can be expressed as a simple time-of-flight equation:t=t_Tx(r,q)+t_Rx(r,p).

The measurement hypersurface corresponding to the above equation is then:

$\left(t,p\right)\in \hat{\Gamma }\left(r,q\right)\iff \{\begin{array}{c}p\in {\prod }^i\\ t=t_{\mathrm{Tx}}\left(r,q\right)+t_{\mathrm{Rx}}\left(r,p\right)\end{array}.$

The parametric formulation of the adjoint operator of the measurement model is next discretised. The same process is applied as for the discretisation of the measurement model. In this case, the discretisation process is more direct since it is directly applied onto p. The discretised set of parameters is Π_dⁱ. We thus have the following equation:{circumflex over (γ)}(r)=Σ_j=1^NeltΣ_i=1^NelrΣ_qk∈XdjΣ_pn∈Πdi(m(ξⁱ)*u)(t_Tx(r,q^k)+t_Rx(r,pⁿ))o_d(r,pⁿ,q^k)w(q^k)z(pⁿ),where w(q^k) and z(pⁿ) are integration weights, m(ξⁱ)=(m(ξⁱ, t^l))_tl∈Td and u∈^Nt is the matched filter of the pulse shape evaluated at each time sample. It can be noticed that the above equation comprises a set of convolved measurement samples ((m(ξⁱ)*u)(t_Tx(r, q^k)+t_Rx(r, pⁿ)))_pn∈Πdi, referred to as hypersurface measurement samples or values, which line on the measurement hypersurface.

An interpolation of the hypersurface measurement samples to substantially coincide with the measurement grid points is carried out next. The values of the hypersurface measurement samples are usually unknown since they do not lie on the measurement grid points. Thus, they need to be approximated from the measurement samples defined on the measurement grid by means of an interpolation kernel. Let us introduce the interpolation kernel ψ: → such that:(m(ξⁱ)*u)(t_Tx(r,q^k)+t_Rx(r,pⁿ))≈Σ_l=1^Nt(m(ξⁱ)*u)(t^l)ψ(t^l−t_Tx(r,q^k)+t_Rx(r,pⁿ)).

The above equation relates the hypersurface measurement samples (m(ξⁱ)*u)(t_Tx(r, q^k)+t_Rx(r, pⁿ)) to the convolved measurement samples (m(ξⁱ)*u)(t^l).

After this, the discretised and interpolated hypersurface measurement samples are summed. The interpolation equation is inserted into the discretised adjoint operator of the measurement model and the following summation is performed:{circumflex over (γ)}(r)=Σ_j=1^NeltΣ_i=1^NelrΣ_qk∈XdjΣ_pn∈Πdi(Σ_l=1^Nt(m(ξⁱ)*u)(t^lψ(t^l−t_Tx(r,q^k)−t_Rx(r,pⁿ)))o_d(r,pⁿ,q^k)w(q^k)z(pⁿ).

The output of the summation is the estimate of the reflectivity value at the considered reflectivity grid point.

The image reconstruction procedure is explained next. The inverse problem is defined first. The measurement values m can be used to generate an estimate of the reflectivity values {circumflex over (γ)} using the adjoint operator of the measurement model as described above. Similarly, an estimate of the reflectivity values {circumflex over (γ)} can be used to generate an estimate of the measurement values {circumflex over (m)}, which are not equal to m, using the measurement model described earlier. Thus, we now know the relationships between the estimates of the measurement and reflectivity values.

The image reconstruction method aims to retrieve an estimate {circumflex over (γ)} of the reflectivity values evaluated at the reflectivity grid points, given the set of measurements values m. The relationship between the measurement values and the unknown reflectivity values is defined by the following linear inverse problem:

$m=\left(\begin{array}{c}\left\{\gamma \right\}\left(p^1,t^1\right)\\ ⋮\\ \left\{\gamma \right\}\left(p^{N_{\mathrm{el}}},t^{N_t}\right)\end{array}\right),$ which corresponds to the measurement model estimated at each point of the measurement grid.

Moreover, since the operator {γ}(pⁿ, t^l) is linear with respect to γ (can be deduced from the formulae described when evaluating the measurement model at a given grid point), there exists a matrix H∈^NelNt×Nγ (associated with the linear operator ) such thatm=Hγ, which defines the inverse problem.

The optimisation problem associated with the inverse problem is explained next. In order to solve the above inverse problem, it is expressed as the following optimisation problem:

$\underset{\gamma }{\min }\left(\gamma ,m\right)=\left(H\gamma ,m\right)+\lambda \left(\gamma \right),$ where (γ, m) denotes the objective function involved in the optimisation problem, (Hγ, m) is a lower semicontinuous functional accounting for the data discrepancy term, measuring the distance between the estimate of the measurement values Hγ and the measurement values m, (γ) is a lower semicontinuous functional describing an optional prior term, which accounts for additional information such as a specific statistical behaviour, on the reflectivity values and λ>0 is a regularisation parameter.

Some examples of the functionals used for the data discrepancy term are listed below:

• • Euclidean distance: (Hγ, m)=∥Hγ−m∥₂²;
  • Indicator function on the -ball of radius ϵ:

$i_{ϵ}\left(\gamma \right)=\{\begin{array}{c}0\mathrm{if}{H\gamma -m}_2<ϵ\\ +\infty \mathrm{otherwise}\end{array};$

• • -norm: (Hγ, m)=∥Hγ−m∥₁.

Some examples of the functionals used for the prior term are listed below:

• • -norm to the power of p: (γ)=∥γ∥_p^p,
    • when p<1, it is a measure of sparsity of γ.
    • when p∈[1,2], it is a measure of how well the data fit to the generalised Gaussian distribution (GGD).
  • -norm to the power of p in a given model Ψ: (γ)=∥Ψγ∥_p^p
    • Ψ can be a general transform such as the wavelet transform or the Fourier transform.
    • Ψ can be a learned dictionary.

It is next explained how the optimisation problem can be solved with one image prior. Depending on the properties of and R, different types of algorithms may be used.

A case where and are differentiable is explained next. In this case, the inverse problem is solved by finding the roots of the derivative of the objective function (γ, m) with respect to γ, which can be written as

$\frac{\partial }{\partial \gamma }+\lambda \frac{\partial \left(\gamma \right)}{\partial \gamma }=0.$

It can be seen that the problem involves the computation of H in the calculation of the derivative of (Hγ, m). Depending on the function (Hγ, m), the calculation of the derivative may also involve the computation of the adjoint operator of the measurement model H. One famous example of such a problem is the Tikhonov regularisation, where (γ)=∥γ∥₂² and (Hγ, m)=∥Hγ−m∥₂². In this case, we have the following solution: {circumflex over (γ)}=(I+H^†H)⁻¹H^†m, where I∈^Nγ×Nγ is the identity matrix, the N_γ×N_γ square matrix with ones on the main diagonal and zeros elsewhere.

A case where only is differentiable and is convex is explained next. One popular group of methods used to solve such problems are called projected gradient methods. Projected gradient methods exploit the fact that the solution of the optimisation problem satisfies a fixed-point equation. One popular algorithm is forward-backward splitting which involves the following system of fixed-point equations, limited to one equation in this specific example:

$\hat{\gamma }={\mathrm{prox}}_{\mathrm{\tau ℛ}}\left(\hat{\gamma }-\tau \frac{\partial \left(H\gamma ,m\right)}{\partial \gamma }\left(\hat{\gamma },m\right)\right),$ where

${\mathrm{prox}}_{\mathrm{\tau ℛ}}\left(x\right)=\underset{y}{\mathrm{argmin}}\frac{1}{2}{y-x}_2^2+\mathrm{\tau ℛ}\left(x\right)$ denotes the proximity operators associated with , the proximity operator being introduced by P. Combettes and J.-C. Pesquet, “Proximal Splitting Methods in Signal Processing”, Fixed-Point Algorithms for Inverse Problems in Science and Engineering, p. 185-212, 2011. The above equation is denoted as a fixed-point equation since {circumflex over (γ)} appears on both sides of the equation. Thus, the problem is solved by carrying out the following iterations:

$={\mathrm{prox}}_{\mathrm{\tau ℛ}}\left(-\tau \frac{\partial \left(H\gamma ,m\right)}{\partial \gamma }\left(,m\right)\right)$ until a convergence criterion is reached.

Another popular algorithm is the fast iterative shrinkage-thresholding algorithm (FISTA), introduced by A. Beck and M. Teboulle, “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems”, SIAM Journal on Imaging Science, vol. 2, p. 183-202, 2009, which is a variant of the method described above.

One example of the forward-backward splitting is when (γ)=∥γ∥₁ and (Hγ, m)=∥Hγ−m∥₂². In this case, the following iteration is performed:=soft(−2τⁿH^†(H−m),τⁿ),where τⁿ∈ is a hyper-parameter.

Indeed, the proximity operator associated with the -norm is defined as prox_τ∥·∥1 (x)=soft(x, τ), where the soft-thresholding operator is defined element-wise as: sign(x_i) max(|x_i|−τ, 0), where sign(x_i) is equal to 1 when x_i>0, 0 when x_i=0 and −1 otherwise.

A case where and are not differentiable is explained next. In this case, several methods can be used, such as alternating direction method of multipliers (ADMM) introduced by S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers”, Foundations and Trends in Machine Learning, vol. 3, p. 1-122, 2011, or primal-dual forward-backward (PDFB) algorithm introduced by P. Combettes, L. Condat, J.-C. Pesquet and BC Vu, “A forward-backward view of some primal-dual optimization methods in image recovery”, Proceedings of the 2014 IEEE International Conference on Image Processing, p. 4141-4145, 2014.

It is next explained how the optimisation problem with M image priors can be solved. In this case it can be considered that the prior term can be written as:(γ)=Σ_k=1^M_k(γ).

In this case, several algorithms can be used such as: primal-dual forward-backward, parallel proximal algorithm (PPXA) introduced by N. Pustelnik, C. Chaux and J.-C. Pesquet, “Parallel Proximal Algorithm for Image Restoration Using Hybrid Regularization”, IEEE Transactions on Image Processing, vol. 20, p. 2450-2462, 2011, or simultaneous direction method of multipliers (SDMM) introduced by S. Setzer, G. Steidl and T. Teuber, “Deblurring Poissonian Images by Split Bregman Techniques”, Journal of Visual Communication and Image Representation, vol. 21, p. 193-199, 2010. The common properties of these algorithms are that they require the computation of the proximity operator of (Hγ, m) (or the derivative of (Hγ, m)) and of the M image prior terms. The computation of the derivative or the proximity operator of (Hγ, m) may involve the estimation of the measurement model and the adjoint operator of the measurement depending on the function (Hγ, m).

The above described image reconstruction method can be summarised by the flow charts of FIGS. 5 to 9. The flow chart of FIG. 5 gives an overview of the method. In step 21, a set of pulse waves is transmitted by a set of transmitting elements. In this example, the transmitting elements have a finite size and arbitrary shape. The set of transmitted pulses may comprise only one pulse or more than one pulses transmitted in a sequence. In step 23, a set of reflected echo waveforms is received by a set of sensors, which in this example have a finite size and arbitrary shape. In step 25, the inverse problem relating the set of measurement values to the set of reflectivity values is defined. In step 27, the solution of the inverse problem is expressed as an optimisation problem involving an objective function comprising a data discrepancy term and one or more image prior terms. In step 29, the optimisation problem is solved by using projections to obtain an estimate of the reflectivity values.

The flow chart of FIG. 6 outlines the proposed method in more detail in the case where the functional is differentiable. The first four steps, ie steps 31, 33, 35 and 37 are the same as the first four steps in the flow chart of FIG. 6. In step 39, a first reflectivity value estimate =H^†m is computed by estimating the adjoint operator of the measurement model on the measurement values m as an initialisation. In step 41, the variable vⁿ= is computed or estimated according to the flow chart of FIG. 8. In step 43, the derivative

$\frac{\partial }{\partial }\left(v^n,m\right)$ is computed. In certain implementation examples, in step 45, the adjoint operator is used to compute the derivative of the data discrepancy term. For instance, when (Hγ, m)=∥Hγ−m∥₂², the derivative of (Hγ, m) is expressed as:

$\frac{\partial }{\partial }\left(v^n,m\right)=H^{†}\left(v^n-m\right),$ which involves the computation of the adjoint operator of the measurement model. In step 47, the proximity operators of the image prior terms are computed. In step 49, is updated. In step 50, it is determined whether or not a convergence criterion is fulfilled. In the affirmative, the process comes to an end. If the criterion is not fulfilled, then the process continues in step 41.

The flow chart of FIG. 7 outlines the proposed method in more detail in the case where the data discrepancy term is the Euclidean distance and where the inverse problem contains one image prior term. The first five steps, ie steps 51, 53, 55, 57, 59 are substantially the same as the first five steps in the flow chart of FIG. 6 with the difference that now the inverse problem comprises one image prior term. In step 61, the variable vⁿ= is computed for each measurement grid point according to the flow chart of FIG. 8. In step 63, the residuals rⁿ=vⁿ−m are computed. In step 65, the adjoint operator of the measurement model on the residuals H^†rⁿ is estimated according to the flow chart of FIG. 10. In step 67, the value −τⁿH^†rⁿ is computed. In step 69, the projection on the image prior term =(−τⁿH^†rⁿ) is computed. In step 70, it is determined whether or not a convergence criterion is fulfilled. In the affirmative, the process comes to an end. If the criterion is not fulfilled, then the process continues in step 61.

The flow chart of FIG. 8 explains an example process of estimating or processing the measurement model. In step 71, a set of reflectivity parametric equations are generated for each measurement grid point (defining a finite size receiving sensor of arbitrary shape), for each point source location in the finite size receiver of arbitrary shape corresponding to the measurement grid point, for each finite size transmitter of arbitrary shape and for each point source location in the finite size transmitter of arbitrary shape. In other words, a set of reflectivity parametric equations is generated for each transmitter-receiver wave path between every pair of the point sources in each transmitter and each receiver. The reflectivity parametric equations define a set of reflectivity hypersurfaces. In step 73, the reflectivity hypersurface is discretised to obtain a set of hypersurface reflectivity samples. In step 75, the hypersurface reflectivity samples are interpolated on the reflectivity grid. In step 77, the interpolated and discretised hypersurface reflectivity samples are summed. In step 79, the value obtained in step 77 is added to the measurement value estimate at the point of the measurement grid. In step 80, the obtained accumulated measurement value estimate is convolved with the pulse shape.

The flow chart of FIG. 9 explains an example process of estimating or processing the adjoint operator of the measurement model. In step 81, the measurement value estimates are convolved with the matched filter of the pulse shape. In step 83, a set of measurement parametric equations are generated for each point of the reflectivity grid, for each finite size transmitter of arbitrary shape, for each finite size receiver of arbitrary shape, for each point source location in the finite size transmitter of arbitrary shape and for each point source location in the finite size receiver of arbitrary shape. In other words, a set of measurement parametric equations is generated for each transmitter-receiver path between every pair of the point sources in each transmitter and each receiver. The measurement parametric equations define a set of measurement hypersurfaces. In step 85, the measurement hypersurface is discretised to obtain a set of hypersurface measurement samples. In step 87, the hypersurface measurement samples are interpolated on the measurement grid. In step 89, the interpolated and discretised hypersurface measurement samples are summed. In step 90, the value obtained in step 89 is added to the reflectivity value estimate at the point of the reflectivity grid.

The transmitted pulse waves can be generated by at least one electromechanical transduction device excited by an electrical signal. The measurement values can be obtained from an electrical signal generated by at least one reciprocal electromechanical transduction device. The electromechanical transduction devices are spatially arranged into: one or more linear array probes with a multi-element structure with a substantially one-wavelength-pitch, and aligned along a straight, convex or concave line; one or more phased array probes with a multi-element structure with a substantially half-wavelength-pitch, and aligned along a straight, convex or concave line; one or more matrix array probes with a multi-element structure with a substantially one-wavelength-pitch, and aligned on a planar surface; or one or more matrix array probes with a multi-element structure with a substantially half-wavelength-pitch, and aligned on a planar surface.

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 transmit acoustic waves are replaced with electromagnetic-radiation pulses. In this particular case, the whole formalism introduced in the present description can be simplified by setting all acoustic transmission propagation delays to zero, because electromagnetic radiation propagates at a speed close to the speed of light, ie typically five orders of magnitude faster than the speed of sound. For this reason, the transmission propagation delays can thus be neglected when accounting the acoustic propagation delays on reception. Another modification considered in the present invention is in the use of transducer technologies different from piezoelectric materials, such as capacitive machined ultrasound transducers (CMUTs). Yet another modification is to apply the invention in a time-multiplexed scheme, whereby subsets of transmitting and/or receiving transducer elements are successively addressed in sequence, and the resulting signals are stored for subsequent recombination. All the above applications can also be translated to microbubble based contrast specific nonlinear imaging techniques known in the art, such as pulse inversion, amplitude modulation, amplitude modulation pulse inversion, harmonic imaging etc.

The present invention also proposes a computer program product comprising instructions stored on a non-transitory medium for implementing the steps of the method as explained above when loaded and run on computing means of an electronic data processing device, such as the imaging apparatus 1.

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.

Claims

1. A pulse wave image reconstruction method comprising: transmitting a set of pulse waves to an object to be at least partially reconstructed by a set of transmitting elements;receiving a set of reflected echo waveforms, by a set of sensors, from the object characterised by a set of reflectivity values defined on a reflectivity grid comprising reflectivity grid points, the set of echo waveforms defining a set of measurement values defined on a measurement grid comprising measurement grid points;defining an inverse problem for obtaining a set of reflectivity value estimates from the set of measurement values by a measurement model;expressing the inverse problem as an optimisation problem involving an objective function comprising a data discrepancy term;expressing the optimisation problem as a system of equations involving contributions of the data discrepancy term;obtaining the set of reflectivity value estimates by performing iterations in the system of equations until a given convergence threshold has been reached, at least some of the iterations involving at least a measurement value estimate for each measurement grid point of a set of measurement grid points and a reflectivity value estimate for each reflectivity grid point of a set of reflectivity grid points;wherein obtaining the contributions in the system of equations comprises:obtaining the measurement value estimate for each measurement grid point of the set of measurement grid points from the measurement model, comprising: generating a set of reflectivity parametric equations per transmitting element-sensor path defining a set of reflectivity hypersurfaces comprising hypersurface reflectivity values;integrating the hypersurface reflectivity values over each reflectivity hypersurface of the set of reflectivity hypersurfaces to obtain the measurement value estimate for the respective measurement grid point;obtaining the reflectivity value estimate for each reflectivity grid point of the set of reflectivity grid points from an adjoint operator of the measurement model, comprising: generating a set of measurement parametric equations per transmitting element-sensor path defining a set of measurement hypersurfaces comprising hypersurface measurement values;integrating the hypersurface measurement values over each measurement hypersurface of the set of measurement hypersurfaces to obtain the reflectivity value estimate for the respective reflectivity grid point.

2. The method according to claim 1, wherein each transmitting element of the set of transmitting elements is considered as a set of transmitting point sources, while each sensor of the set of sensors is considered as a set of receiving point sources.

3. The method according to claim 1, wherein, before integrating the hypersurface reflectivity values, the method further comprises: discretising each reflectivity hypersurface of the set of reflectivity hypersurfaces to obtain the hypersurface reflectivity values;interpolating the discretised hypersurface reflectivity values to substantially coincide with the reflectivity grid points.

4. The method according to claim 1, wherein, before integrating the hypersurface measurement values, the method further comprises: discretising each measurement hypersurface of the set of measurement hypersurfaces to obtain the hypersurface measurement values;interpolating the discretised hypersurface measurement values to substantially coincide with the measurement grid points.

5. The method according to claim 1, wherein the data discrepancy term is expressed as a positive functional evaluating a distance between the measurement value and the measurement value estimate for each measurement grid point of the set of measurement grid points.

6. The method according to claim 5, wherein the positive functional is a positive differentiable functional.

7. The method according to claim 1, wherein the objective function is expressed as a combination of the data discrepancy term and at least one image prior term, the image prior term comprising prior information on the reflectivity values.

8. The method according to claim 7, wherein the at least one image prior term is expressed as the 1p-norm to the power of p of the representation of the set of reflectivity value estimates in a given model, the 1p-norm being defined for a vector a as α∈N as

9. The method according to claim 1, wherein the pulse waves are pulse acoustic waves or electromagnetic-radiation waves.

10. The method according to claim 1, wherein the pulse waves comprise a set of steered plane waves or a set of diverging waves.

11. The method according to claim 1, wherein the pulse waves are generated by at least one electromechanical transduction device excited by an electrical signal, and the measurement values are obtained from an electrical signal generated by at least one reciprocal electromechanical transduction device.

12. The method according to claim 11, wherein the electromechanical transduction devices are spatially arranged into a linear array probe with a multi-element structure with a substantially one-wavelength-pitch, and aligned along a straight, convex or concave line;a phased array probe with a multi-element structure with a substantially half-wavelength-pitch, and aligned along a straight, convex or concave line;a matrix array probe with a multi-element structure with a substantially one-wavelength-pitch, and aligned on a planar surface; ora matrix array probe with a multi-element structure with a substantially half-wavelength-pitch, and aligned on a planar surface.

13. The method according to claim 1, wherein the reflectivity and measurement hypersurfaces have a conical shape when the set of reflectivity values lie in a 2D space, and a quadric shape when the set of reflectivity values lie in a 3D space.

14. The method according to claim 1, wherein the integration of the hypersurface measurement or reflectivity values is achieved by numerical integration, and wherein the method further comprises: adding the obtained measurement value estimate to the previous measurement value estimate at the respective measurement grid point to obtain an accumulated measurement value estimate; and convolving the accumulated measurement value estimate with the pulse shape of the transmitted pulse wave.

15. An imaging apparatus for reconstructing an image of an object, the apparatus being configured to perform operations comprising: transmit a set of pulse waves to an object to be at least partially reconstructed by a set of transmitting elements;receive a set of reflected echo waveforms, by a set of sensors, from the object characterised by a set of reflectivity values defined on a reflectivity grid comprising reflectivity grid points, the set of echo waveforms defining a set of measurement values defined on a measurement grid comprising measurement grid points;define an inverse problem for obtaining a set of reflectivity value estimates from the set of measurement values by a measurement model;express the inverse problem as an optimisation problem involving an objective function comprising a data discrepancy term;express the optimisation problem as a system of equations involving contributions of the data discrepancy term;obtain the set of reflectivity value estimates by performing iterations in the system of equations until a given convergence threshold has been reached, at least some of the iterations involving at least a measurement value estimate for each measurement grid point of a set of measurement grid points and a reflectivity value estimate for each reflectivity grid point of a set of reflectivity grid points;wherein for obtaining the contributions in the system of equations, the apparatus is configured to perform operations comprising:obtain the measurement value estimate for each measurement grid point of the set of measurement grid points from the measurement model, the measurement value estimate being configured to be obtained by: generating a set of reflectivity parametric equations per transmitting element-sensor path defining a set of reflectivity hypersurfaces comprising hypersurface reflectivity values;integrating the hypersurface reflectivity values over each reflectivity hypersurface of the set of reflectivity hypersurfaces to obtain the measurement value estimate for the respective measurement grid point;obtain the reflectivity value estimate for each reflectivity grid point of the set of reflectivity grid points from an adjoint operator of the measurement model, the reflectivity value estimate being configured to be obtained by:generating a set of measurement parametric equations per transmitting element-sensor path defining a set of measurement hypersurfaces comprising hypersurface measurement values;integrating the hypersurface measurement values over each measurement hypersurface of the set of measurement hypersurfaces to obtain the reflectivity value estimate for the respective reflectivity grid point.