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 realtime 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:

• Figure 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;

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

• Figure 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;

• Figure 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; · Figure 5 is a flow chart summarising the image reconstruction method according to an example of the present invention;

• Figure 6 is a flow chart illustrating the image reconstruction method in more detail according to an example of the present invention; • Figure 7 is a flow chart illustrating the image reconstruction method in more detail according to another example of the present invention;

• Figure 8 is a flow chart illustrating a process of implementing an image measurement model according to an example of the present invention; and · Figure 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 = Ηγ + 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 postprocessing and display modules, as described in Figure 1 .

Figure 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 1 D 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 of the image under scrutiny. The estimate can be written as:

where / 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, E) has an explicit formulation (eg backpropagation methods using the adjoint operator of H).

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

where J(y, m) denotes an objective function which comprises a functional T, also referred to as a data discrepancy term, measuring the error between 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 postprocessing and display unit or module 1 1 . 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) y(r) with r = [x, y, z] ^T E Ω are detected by sensors 13. Figure 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.

• denotes the -norm of a vector where R

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

is defined as

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).

• y(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. · α(r, p) defines the value of the attenuation between a point source located at p and a point located at r. • is tne 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: where refers to the number of

transmitting elements, also referred to as transmitters. where refers to the number of

receive sensors, also referred to as receivers.

• A set of point source locations inside a transmit element:

· A set of point source locations inside a receive sensor:

• A set of time samples:

Reflectivity grid related notations:

• A set of locations of reflectivity values:

where Ν _γ 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:

The above equation can be expressed as shown below, and defines the measurement model as an integral along the hypersurface Γ(ρ, q, t):

where

In the above equations, the hypersurface, or more specifically a reflectivity hypersurface, is Γ(ρ, q, t), defined by the equation g(r, p, q, t) = 0, where the variable of interest is r. 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. Figure 3 schematically illustrates the measurement configuration in continuous domain showing also the reflectivity hypersurface.

The reflectivity hypersurface r(p, q, t) can now be parametrised in order to reduce the integration dimension. The reflectivity hypersurface Γ(ρ, 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:

The equation on the right of the equivalence defines the reflectivity parametric equations of r since r now depends on a parameter α E A, where

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: where

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* . The measurement model is defined as

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:

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 It can be for example the coordinate of

the point in the centre of Π ¹ .

It can be written: 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:

The corresponding parametric formulation is given by:

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:

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:

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 as

which is only defined for discrete coordinates In the remainder of the present description, the following holds:

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 of the parameter α, which defines the following set of parameters: A _d =

It can be seen that the selection of the discrete values of the

parameter α leads to a discretisation of the hypersurface which will now be

described by only N _a points Figure 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:

where is the pulse shape evaluated at each time sample of the measurement

where 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 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. Thevenaz, 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:

The above interpolation equation relates the hypersurface reflectivity values

(which are typically not on the reflectivity grid points) to the

reflectivity samples

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:

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 such that

then the adjoint operator is

defined as:

Let us consider the following inner product:

wher s the matched filter of the pulse echo waveform. We now have:

where

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:

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:

The measurement hypersurface corresponding to the above equation is then:

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 We thus have the

following equation:

where are integration weights,

N ^t 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 referred to as hypersurface

measurement samples or values, which lie 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 kerne such that:

The above equation relates the hypersurface measurement samples to the convolved measurement samples

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:

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 using the adjoint operator of the measurement model as described above. Similarly, an estimate of the reflectivity values can be used to generate an

estimate of the measurement values 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 γ 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:

which corresponds to the measurement model estimated at each point of the measurement grid.

Moreover, since the operator s linear with respect to y (can

be deduced from the formulae described when evaluating the measurement model at a given grid point), there exists a matrix (associated with the linear

operator such that

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:

where denotes the objective function involved in the optimisation problem, is a lower semicontinuous functional accounting for the data discrepancy

term, measuring the distance between the estimate of the measurement values Ηγ and the measurement values 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:

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

• when p < 1, it is a measure of sparsity of γ.

• when 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 mode

· Ψ 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 different types of algorithms may

be used.

A case where are differentiate is explained next. In this case, the

inverse problem is solved by finding the roots of the derivative of the objective function with respect to γ, which can be written as

It can be seen that the problem involves the computation of H in the calculation of the derivative of Depending on the function 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 In this case, we

have the following solution where is the identity

matrix, the Ν _γ χΝ _γ square matrix with ones on the main diagonal and zeros elsewhere.

A case where only T 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:

where denotes the proximity operators

associated with Jl, 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 appears on both sides of the equation.

Thus, the problem is solved by carrying out the following iterations:

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 ^an d

In this case, the following iteration is performed:

where s a hyper-parameter.

Indeed, the proximity operator associated with the -norm is defined as

where the soft-thresholding operator is defined element-

is equal to 1 when x _t > 0, 0 when

A case where T and JZ 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, 201 1 , 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:

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, 201 1 , 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 (or the derivative of and of the M image

prior terms. The computation of the derivative or the proximity operator of

may involve the estimation of the measurement model and the adjoint operator of the measurement depending on the function

The above described image reconstruction method can be summarised by the flow charts of Figures 5 to 9. The flow chart of Figure 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 Figure 6 outlines the proposed method in more detail in the case where the functional T is differentiate. The first four steps, ie steps 31 , 33, 35 and 37 are the same as the first four steps in the flow chart of Figure 6. In step 39, a first reflectivity value estimate is computed by estimating the adjoint

operator of the measurement model on the measurement values m as an initialisation. In step 41 , the variable is computed or estimated according to the flow chart

of Figure 8. In step 43, the derivative ' ^s 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

the derivative of is expressed as: 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 Figure 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 Figure 6 with the difference that now the inverse problem comprises one image prior term. In step 61 , the variable is computed for each measurement grid point according to

the flow chart of Figure 8. In step 63, the residuals r ⁿ = v ⁿ —m are computed. In step 65, the adjoint operator of the measurement model on the residuals is estimated according to the flow chart of Figure 10. In step 67, the value is computed.

In step 69, the projection on the image prior term 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 Figure 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 Figure 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.