SIGNAL SAMPLING WITH JOINT TRAINING OF LEARNABLE PRIORS FOR SAMPLING OPERATOR AND DECODER

Abstract

A method of sampling and decoding of a signal of interest x comprising, at a training stage: acquiring a set of training signals {x.sub.i}.sub.i=1.sup.M, providing a sampling operator P.sub. and a decoder g.sub.g(.), training operator P.sub. on signals {x.sub.i}.sub.i=1.sup.M to obtain a learned sampling operator P.sub.{circumflex over ()}; and, at a sampling stage: applying P.sub.{circumflex over ()} in a transform domain If to signal x, resulting in observation signal y; applying the decoder g.sub.g(.) to y, to produce an estimate {circumflex over (x)} of signal x to decode and/or, decide about, the signal. Decoder g.sub.g(.) is trained jointly with P.sub. on signals {x.sub.i}.sub.i=1.sup.M, to obtain a learned decoder g.sub.g, by jointly determining, during a cost minimization step, sampling parameters and decoding parameters .sub.gaccording to a cost function, and wherein the step of applying the decoder g.sub.g(.) uses decoding parameters .sub.g, such that estimate {circumflex over (x)} is produced by the learned decoder g.sub.{circumflex over ()}g.

Claims

1. A method of sampling and decoding of a signal of interest x comprising: at a training stage, the steps of: acquiring a set of training signals {x.sub.i}.sub.i=1.sup.M, providing a sampling operator P.sub. and a decoder g.sub.g(.), training said sampling operator P.sub. on said set of training signals {x.sub.i}.sub.i=1.sup.M to obtain a learned sampling operator P.sub.{circumflex over ()}, training said decoder g.sub.g(.) jointly with said sampling operator P.sub. on said set of training signals {x.sub.i}.sub.i=1.sup.M, to obtain a learned decoder g.sub.{circumflex over ()}g, by jointly determining, during a cost minimization step of the training stage, a set of sampling parameters and a set of decoding parameters .sub.g according to a cost function, and at a sampling stage, the steps of: applying during a sampling step the learned sampling operator P.sub.{circumflex over ()} in a transform domain to the signal of interest x, resulting in observation signal y, applying during a decoding step the decoder g.sub.g(.) to the observation signal y, to produce an estimate {circumflex over (x)} of the signal of interest x in order to decode the signal and/or to make a decision about the signal, wherein applying the decoder g.sub.g(.) includes using said set of decoding parameters .sub.g, such that the estimate {circumflex over (x)} of the signal of interest x is produced by the learned decoder g.sub.{circumflex over ()}g.

2. The method according to claim 1, wherein the cost minimization step for determining the set of sampling and decoding parameters .sub.g during the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) is implemented by an optimization procedure iterated in alternative directions until convergence.

3. The method according to claim 2, wherein said optimization procedure iterated in alternative directions until convergence is based on: an optimization of the estimate {circumflex over (x)} of the signal of interest x and of the decoding parameters .sub.g, whilst keeping the set of sampling parameters fixed, of the form
{circumflex over (x)}.sub.=argmin.sub.xyP.sub.x.sup.2.sub.2+.sub.x.sub.x(x), where .sub.x(x) is a signal regularization prior, .sub.x is a regularization parameter and ..sub.2 denotes the l.sub.2-norm, and on an optimization of the set of sampling parameters , whilst keeping the decoding parameters .sub.gfixed, of the form
{circumflex over ()}=argmin.sub..sub.i=1.sup.Mx.sub.i{circumflex over (x)}.sub.i.sup.2.sub.2+.sub..sub.(), where .sub.(.) is a prior on desirable properties of the set of sampling parameters .

4. The method according to claim 2, wherein said optimization procedure iterated in alternative directions until convergence is based on an optimization of the general form:
({circumflex over ()}, {circumflex over ()}.sub.g)=argmin.sub.,.sub.g.sub.p=1.sup.P.sub.i=1.sup.Mx.sub.ig.sub.g(P.sub.(x.sub.i+z.sub.p)).sup.2.sub.2+.sub.1.sub.()+.sub.2.sub.g(.sub.g), where the decoder g.sub.g(.) is factorized as a deep structure of a form g.sub.g(x)=.sub.D.sup.g(W.sub.D.sup.g . . . .sub.1.sup.g(W.sub.1.sup.g(x))), with .sub.g=(W.sub.1.sup.g . . . , W.sub.D.sup.g) being parameters of a deep network, .sub.k.sup.g(.) standing for point-wise non-linearities k=1, . . . , D with D denoting the number of layers, .sub.() denoting constraints on the geometry and properties of the sampling operator P.sub., .sub.g(.sub.g) defining constraints on the decoder parameters, .sub.1 and .sub.2 being regularization parameters and ..sub.2 denoting the l.sub.2-norm.

5. The method according to claim 4, wherein said optimization procedure iterated in alternative directions until convergence is based on an optimization of the decoding parameters .sub.g, whilst keeping the set of sampling parameters fixed, of the form
{circumflex over ()}.sub.g=argmin.sub..sub.g.sub.i=1.sup.Mx.sub.ig.sub.g(P.sub.(x.sub.i)).sup.2.sub.2+.sub.1.sub.g(.sub.g), and on an optimization of the set of sampling parameters , whilst keeping the decoding parameters .sub.g fixed, of the form
{circumflex over ()}=argmin.sub..sub.i=1.sup.Mx.sub.i{circumflex over (x)}.sub.i.sup.2.sub.2+.sub.1.sub.().

6. The method according to claim 1, wherein the cost function used in the cost minimization step for determining the set of sampling and decoding parameters .sub.g during the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) is of the general form ({circumflex over ()},{circumflex over ()}.sub.g)=argmin.sub.,.sub.g.sub.p=1.sup.P.sub.i=1.sup.M.sub.m=1.sup.CL(l(m),g.sub.g(P.sub.(x.sub.i+z.sub.p)))+.sub.1.sub.()+.sub.2.sub.g(.sub.g), wherein the decoder g.sub.g(.) is a classifier g.sub.g:.sup.n.fwdarw.{1, . . . , C}, .sub.() denotes constraints on the geometry and properties of the sampling operator .sub., .sub.g(.sub.g) denotes constraints on the decoder parameters, .sub.1 and .sub.2 being regularization parameters and l(m) is an encoded class label for the class g.sub.g:.sup.n.fwdarw.{1, . . . , C} with L(.) denoting the cost function.

7. The method according to claim 1, wherein both the sampling operator P.sub. and the decoder g.sub.g(.) are realized on multiple stages.

8. The method according to claim 1, wherein the sampling operator P.sub. is repartitioned into several subbands which are optimized to the properties of the set of training signals {x.sub.i}.sub.i=1.sup.M and of the decoder g.sub.g(.) according to a chosen cost function and/or into several transform domains .

9. The method according to claim 1, wherein the sampling operator P.sub. is specifically adapted for a transform domain representing a color space transformation by performing the sampling in each color component of said color space transform domain .

10. The method according to claim 1, wherein the sampling operator P.sub. has a reduced sampling rate such as to ensure privacy-preserving sampling of the signal of interest x in remote or on-device authentication and identification applications.

11. The method according to claim 1, wherein the sampling operator P.sub. is designed such as to depend on both the set of training signals {x.sub.i}.sub.i=1.sup.M and a secret key and/or a product ID to ensure a secret subset sampling in order to prevent counterfeiting and adversarial attacks.

12. The method according to claim 1, wherein the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) is performed based on a class of training signals corresponding to genuine groups of objects, products or persons and non-genuine ones, optimization of a cost function ensuring an optimal distinguishability of genuine from non-genuine classes.

13. The method according to claim 1, wherein the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) is performed based on a set of training signals {x.sub.i}.sub.i.sub.i=1.sup.M representing the same objects under different viewing and imaging conditions.

14. The method according to claim 1, wherein the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) uses a model of the signal of interest x which is based on an auto-encoder comprising a model encoder and decoder pair and implemented by taking into account compressive and/or robustification constraints in the space of latent representation.

15. The method according to claim 1, wherein the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) uses a model of the signal of interest x which is based on a single learnable orthonormal or overcomplete representation or a deep learning representation which either is also trained based on the set of training signals {x.sub.i}.sub.i=1.sup.M or is trained on the fly.

16. The method according to claim 1, wherein the set of training signals {x.sub.i}.sub.i=1.sup.M comprises both clean and distorted and/or corrupted signals, the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) ensuring selection of the most informative and distinguishable sets of sampling and decoding parameters .sub.g.

17. A computer program means stored in a tangible computer readable medium adapted to implement the method according to claim 1.

18. A device equipped with the computer program means according to claim 17, wherein the device is chosen from a group comprising a mobile phone equipped with a camera, a digital photo apparatus, a digital video camera, a scanning device, a tablet or personal computer, and a server.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0028] The attached drawings exemplarily and schematically illustrate the principles as well as several embodiments of the present invention.

[0029] FIG. 1 presents a generic architecture of training and sampling where the sampling operator and the decoder are jointly trained on training signals and then applied to a test signal.

[0030] FIG. 2 further explains the method of training and sampling of FIG. 1 by showing the details of the main elements of the sampling and decoding pipelines combined with a minimization of a chosen cost function.

[0031] FIG. 3 is a simplified block-diagram illustrating a model training framework based on a set of training signals resulting in the estimation of encoder and decoder model parameters.

[0032] FIG. 4 schematically illustrates an implementation of an iterative decoder in a form of a multi-layer network with fixed blocks of model encoding and decoding as well as common reconstruction blocks A and B chosen to be fixed for illustration purposes.

[0033] FIG. 5 schematically illustrates the difference between signal sampling according to the present invention as compared to two prior art approaches to signal sampling.

[0034] FIG. 6 illustrates a joint training of the sampling operator and the decoder for a generic classification problem.

[0035] FIG. 7 shows a first embodiment of a multi-stage sampling-decoding framework according to the present invention.

[0036] FIG. 8 shows a second embodiment of a multi-stage sampling-decoding framework which is extended as compared to the example of FIG. 7.

[0037] FIG. 9 illustrates application of multi-stage sampling to spatio-temporal signals.

[0038] FIG. 10 presents the results of sampling and reconstruction for one medical (MRI) and two astronomical datasets of images denoted as D1, D2 and D3 for the Euclidean cost function for different sampling rates.

[0039] FIG. 11 illustrates some examples of images tested in order to obtain the results shown in FIG. 10.

DETAILED DESCRIPTION OF THE INVENTION

[0040] In the following, the invention shall be described in detail with reference to the above mentioned figures.

[0041] The present invention relates to a method of sampling of signals based on especially learned sampling operators and decoders as well as to a device implementing such a sampling method. As already mentioned before, the following description, in general, will concentrate on the consideration of the method according to the present invention when used for the above mentioned reconstruction problem, i.e. for image reconstruction and decoding, and will only highlight, respectively exemplify, differences of the method when used for the above mentioned classification problems in the course of the description. Furthermore, as usual in the domain of signal, image and video processing, the term signal is considered throughout the description to cover both analog and digital data in one or more dimensions.

[0042] According to the present invention, the sampling operator is trained in such a way to ensure the most accurate and discriminative estimation of signal model parameters that in turn lead to accurate and efficient decoding for the reconstruction or classification. Accordingly, the method according to the present invention comprises two stages: A first stage includes joint learning of a sampling operator P.sub. and of a generalized decoder g.sub.g(.), which has a set of parameters .sub.g, the first stage being referred to as learning or training stage. A second stage includes testing comprising reconstruction of a signal of interest x from an observed signal y and/or other tasks such as detection, classification, recognition, identification and/or authentication using the previously learned, i.e. trained, sampling operator P.sub.{circumflex over ()} and decoder g.sub.{circumflex over ()}g, the second stage being referred to as sampling or testing stage. In the context of identification and authentication, as part of a general recognition problem, the training stage is also called enrolment and the testing stage is referred to as verification stage.

[0043] FIG. 1 illustrates the principle of these two stages for the general case of application of the method to the above mentioned reconstruction problem. By use of a set of training signals {x.sub.i}.sub.i=1.sup.M(20) acquired during an acquisition step (10), a joint training of a sampling operator P.sub. and of a decoder g.sub.g(.) is performed in training stage (30), resulting into a particular set of sampling parameters {circumflex over ()} used for the learned sampling operator P.sub.(40) and a particular set of decoding parameters {circumflex over ()}.sub.g used for the learned decoder g.sub.{circumflex over ()}g(50). The sampling of unknown signal of interest x(60) is performed in sampling stage (70) using the given, learned sampling operator P.sub.{circumflex over ()}(40) resulting into the observation signal y(80). The following decoding step (90) is based on the given, learned decoder g.sub.{circumflex over ()}g and results into an estimate {circumflex over (x)}(100) of the signal of interest. For reasons of completeness, it shall be noted in this context that, in general, it is possible, depending on the specific application of the method according to the present invention, that the transform domain is chosen to be an identity operator =I, such that the joint training of the sampling operator P.sub. and of the decoder g.sub.g(.) is performed in the spatial domain instead of in some transform spectral domain.

[0044] More specifically, the training procedure generally illustrated in FIG. 1 is based on the training pipeline shown in more detail in FIG. 2. The set of training signals {x.sub.i}.sub.i=1.sup.M is observed via a sequence of transformations using a transform domain (110), the sampling operator P.sub.(120) and an optional quantization Q(.)(130), producing a set of observable signals {y.sub.i}.sub.i=1.sup.M that might be additionally corrupted by noise. These operations are combined into a block by an observation encoder (160). A parametrized observation decoder (170) includes the decoder g.sub.g(140) that produces estimates {{circumflex over (x)}.sub.i}.sub.i=1.sup.M(150) of the set of training signals {x.sub.i}.sub.i=1.sup.M based on the set of observable signals {y.sub.i}.sub.i=1.sup.M . The training is performed, during cost minimization step (180), by minimizing a cost function leading to the best possible reconstruction of the estimates {{circumflex over (x)}.sub.i}.sub.i=1.sup.M closest to {x.sub.i}.sub.i=1.sup.M in some metric or latent space over sets of sampling parameters and of decoding parameters .sub.g. In the simplest case, the Euclidean distance can be used. The results of this optimization are estimates of the learned sampling operator P.sub.{circumflex over ()}(40) and of the learned decoder g.sub.{circumflex over ()}g(50). Analogously, the sampling procedure generally illustrated in FIG. 1 is based on the sampling pipeline shown in more detail in FIG. 2 and combines an observation encoder (160) and decoder (170) with the trained sampling operator P.sub. and decoder g.sub.{circumflex over ()}g used in steps (120) and (140). The result of decoding is an estimate {circumflex over (x)}(100) of the signal of interest x.

[0045] The present invention proposes, in particular in the context of the above mentioned reconstruction problem, two approaches for the joint training stage of the sampling operator and of the decoder as well as for the testing/sampling stage schematically shown in FIG. 2. Like will become clear in the following, the two approaches mainly differ at the training stage, whilst the corresponding sampling stages of the two approaches are almost identical. In the following, these approaches will be referred to as inverse problem solution and regression approaches.

[0046] According to the inverse problem solution approach to reconstruction, the joint training of the sampling operator and of the decoder comprises two steps. During a first step, the set of sampling parameters is fixed and one tries to find the estimate {circumflex over (x)}(100) of the signal of interest and the decoding parameters .sub.g by solving an optimization problem using hand-crafted or learned properties of the signal of interest x from the set of training signals {x.sub.i}.sub.i=1.sup.M. The optimization is described by equation

{circumflex over (x)}.sub.=argmin.sub.xyP.sub.x.sup.2.sub.2+.sub.x.sub.x(x), (4)

where .sub.x(x) is a signal regularization prior linked to the properties and priors of the signals to be sampled and reconstructed, .sub.xis a regularization parameter and ..sub.2denotes the l.sub.2-norm. The regularization prior .sub.x(x) and the regularization parameter .sub.x aren't used in prior art methods and allow to take into account in improved manner the predetermined properties of specific types of signals, i.e. of the kind of images, to be sampled and reconstructed by introducing corresponding priors which are learned during the training stage of the method according to the present invention. In fact, to illustrate utility of these regularization priors and parameters in practical manner, it is clear that medical images, e.g. images of human kidney obtained by MRI, have different, but predetermined properties than images obtained by astronomy observations or else images of a given kind of objects obtained by CMOS/CCD cameras like simple facial portraits or the like. The predetermined properties of the corresponding images are reflected in the method according to the present invention in different, previously trained regularization priors .sub.x(x) and parameters .sub.x for each type of images.

[0047] Different strategies for the design of the signal priors .sub.x(x) will be considered here below. The resulting estimate {circumflex over (x)}(100) of the signal of interest can be expressed by equation

{circumflex over (x)}.sub.=g.sub.{circumflex over ()}g(y,), (5)

where g.sub.{circumflex over ()}g(.) is a decoder obtained from equation (1) and denotes the particular set of decoding parameters. The set of sampling parameters is used in the decoder g.sub.{circumflex over ()}g(.) as an argument to signify the fact that the decoder depends on the set of sampling parameters .

[0048] During a second step, using equation (2) and the set of training signals {x.sub.i}.sub.i=1.sup.M consisting of M training samples, one tries to optimize the set of sampling parameters according to equation

{circumflex over ()}=argmin.sub..sub.i=1.sup.Mx.sub.i{circumflex over (x)}.sub.i.sup.2.sub.2+.sub..sub.(), (6)

where .sub.(.) is a prior on desirable properties of the set of sampling parameters combining the geometry, the number of samples, etc.

[0049] The optimization procedure iterates the first and second steps till convergence. This procedure is similar to an alternative direction minimization.

[0050] After joint training of the sampling operator and of the decoder, the testing/sampling stage of the inverse problem solution approach also has two steps during which the learned sampling operator and decoder are applied to real data, i.e. to a signal, which was not seen before, but which is assumed to follow the same statistics as the training data.

[0051] During a first step, by using a given set of particular sampling parameters {circumflex over ()}, one tries to produce the observation signal y by use of equation

y=P.sub.{circumflex over ()}x. (7)

[0052] During a second step, by using a given observation signal y, the decoder (2) produces an estimate {circumflex over (x)}(100) of the signal of interest using equation

{circumflex over (x)}=g.sub.(y,{circumflex over ()}). (8)

[0053] According to the regression approach to the reconstruction problem, the joint training of the sampling operator and of the decoder may also be formulated as an optimization problem which, in this case, may be described by equation

({circumflex over ()},{circumflex over ()}.sub.g)=argmin.sub.,.sub.g.sub.p=1.sup.P.sub.i=1.sup.Mx.sub.ig.sub.g(P.sub.(x.sub.i+z.sub.p)).sup.2.sub.2+.sub.1.sub.()+.sub.2.sub.g(.sub.g), (9)

where the decoder is assumed to be factorized as a deep structure of a form described by equation g.sub.g(x)=.sub.D.sup.g(W.sub.D.sup.g . . . .sub.1.sup.g(W.sub.1.sup.g(x))), with .sub.g=(W.sub.1.sup.g . . . , W.sub.D.sup.g) being parameters of a deep network, whilst .sub.k.sup.g(.) stands for point-wise non-linearities k=1, . . . , D with D denoting the number of layers; .sub.() denotes the constraints on the geometry and properties of the sampling operator and .sub.g(.sub.g) defines constraints on decoder parameters; .sub.1 and .sub.2 are regularization parameters. Equation (9) also includes a factor for additional robustification to possible noise perturbations via additive component z.sub.p, p=1, . . . , P following some predefined distribution, where P denotes the number of noisy training examples.

[0054] Since the simultaneous minimization is not feasible due to the non-convexity in both and .sub.g, one can proceed by reformulating the problem as an iterative alternative minimization, keeping one of the parameters fixed and optimizing the second one, whilst preferably skipping the noise component for simplicity. Such manner of proceeding also has two steps.

[0055] During a first step, the set of sampling parameters is kept fixed and one tries to find the decoding parameters according to equation

{circumflex over ()}.sub.g=argmin.sub..sub.g.sub.i=1.sup.Mx.sub.ig.sub.g(P.sub.(x.sub.i)).sup.2.sub.2+.sub.1.sub.g(.sub.g). (10)

[0056] During a second step, the decoding parameters {circumflex over ()}.sub.g are kept fixed and one tries to optimize the set of sampling parameters according to equation

{circumflex over ()}=argmin.sub..sub.i=1.sup.Mx.sub.i{circumflex over (x)}.sub.i.sup.2.sub.2+.sub.1.sub.(), (11)

similarly to equation (6) in the inverse problem solution approach. The optimization procedure iterates the first and second steps till convergence.

[0057] The testing/sampling stage of the regression approach is identical to the one of the inverse problem solution approach and is based on above mentioned equations (7) and (8), however with the difference in the structure of the decoder g.sub.g which is represented in the regression approach in a form of a deep network, such as indicated in equation (9).

[0058] For reasons of completeness and of clarity, it is added here that the nomenclature inverse problem solution approach and regression approach pertains to the fact that, on the one hand, in the first named approach the decoder g.sub.{circumflex over ()}g(.) is only implicitly constructed as a solution of equation (4), whilst in the second named approach the decoder g.sub.{circumflex over ()}g(.) is explicitly constructed in equation (9). On the other hand, in the first named approach equation (4) represents a typical signal processing or physics based solution of an inverse problem, this explaining the choice of name for this approach, whilst in the second named approach equation (9) represents a direct joint optimization formulation for both of the sampling operator and of the decoder, also called regression formulation, this being a typical machine learning formulation.

Construction of Priors (x) for the Inverse Problem Solution Approach

[0059] The presence of reliable priors (x) plays a crucial role for the above mentioned reconstruction problem, in particular within the above described inverse problem solution approach, which is mainly due to the, in principle, ill-posed nature of equation (4).

[0060] In general, one can distinguish three cases as far as is concerned the construction of a model for the regularization priors (x) in equation (4), namely a first case where no priors are used, a second case using hand-crafted priors, both of which belong to prior art, and a third case using learnable priors and corresponding to a particularly preferred embodiment of the sampling method according to the present invention. FIG. 5 schematically illustrates the difference between the preferred learning of priors within the signal sampling method according to the present invention as compared to said two prior art approaches to signal sampling. Like symbolically indicated in the left upper part of FIG. 5, the first prior art approach to signal sampling is based on a standard sampling requiring acquisition of all signal samples under the absence of any prior information about the nature of acquired signals. Like symbolically indicated in the right upper part of FIG. 5, the second prior art approach to signal sampling is based on a sampling operator adapted to some fixed linear decoder and chosen metric of reconstruction, thus leading to sampling of the largest magnitude coefficients that are concentrated for natural images in the region of low frequencies. Like symbolically indicated in the lower part of FIG. 5 and though the above mentioned two prior art approaches may in principle also be used within the joint training stage of the method according to the present invention, the method of signal sampling according to the present invention preferably uses learnable encoder and decoder priors.

[0061] To describe these differences in more detail, it is to be noted that the first case corresponds to the method disclosed in the above mentioned article Learning-based subsampling and US 2017/0109650 by V. Cevher et al., such method completely disregarding the above mentioned signal regularization priors, but in fact only solving a least squares problem described by equation

{circumflex over (x)}=argmin.sub.xyP.sub.x.sup.2.sub.2, (12)

leading to a so-called linear decoder described by equation

x=*P.sub..sup.Ty, (13)

that can also be considered as a particular case of equation (4) for .sub.x=0.

[0062] The second case is based on hand-crafted priors assuming that the signal of interest x might be sparsely presented in some orthogonal basis like Fourier, DCT, wavelet, Hadamard, etc. One can distinguish two models here.

[0063] In a first model using a so-called synthesis or sparse approximation model, one can present the signal of interest x by equation

x=u+e.sub.x, (14)

where e.sub.x is a residual error of approximation and u is assumed to be sparse, leading to l.sub.1- or l.sub.0-norm regularizers. The corresponding regularizer can be described by equation

[00001]$\begin{array}{cc}{}_x\left(x\right)=\frac{}{2}.Math.{.Math.x-.Math..Math.u.Math.}_2^2+{}_u.Math.{}_u\left(u\right),& \left(15\right)\end{array}$

where .sub.u(u)is a regularizer which imposes a sparsity on (u) as well as and .sub.u are regularization coefficients. This is described for example by M. Aharon et al. in the article K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation, IEEE Transactions on signal processing, 54, 11, pp. 4311-4322, 2006.

[0064] In contrast, a second model known as the analysis model or also known as sparsifying transform, disclosed for example by S. Ravishankar and Y. Bressler in the article Learning Sparsifying transform, 61, 5-8, pp. 1072-1086, 2013, assumes that

u=Wx+e.sub.u, (16)

where W is the analysis transform that for the orthonormal case is just W=.sup.T and e.sub.u represents the noise of approximation. The corresponding regularizer can be described by equation

[00002]$\begin{array}{cc}{}_x\left(x\right)=\frac{}{2}.Math.{.Math.\mathrm{Wx}-.Math.u.Math.}_2^2+{}_u.Math.{}_u\left(u\right).& \left(17\right)\end{array}$

[0065] The solution to equation (4) with priors according to equations (15) or (17) is generally obtained by iterative methods.

[0066] The third case corresponds to a particularly preferred embodiment of the sampling method according to the present invention and is based on learnable priors that explore the full power of training with the set of training signals {x.sub.i}.sub.i=1.sup.M. In principle, one can reformulate equation (4) jointly with equations (14) and (16) to learn the basis and the analysis transform W as a single overcomplete transform or layer.

[0067] However, in view of the complexity of this problem and the benefits offered by deep architectures, the following will consider an alternative model. Given a set of training signals {x.sub.i}.sub.i=1.sup.M, training of a model prior factorized via a model encoder-decoder pair such as shown in FIG. 3 will be considered as a stand-alone step. For the most general case, it will be assumed that the model encoder .sub.E(.)(300) is trained jointly with the model decoder .sub.D( )(320) according to a predefined cost function (340) to estimate the encoder parameters {circumflex over ()}.sub.E(350) and the decoder parameters {circumflex over ()}.sub.D(360). The parameters of the encoder .sub.E=(W.sub.1.sup.E, . . . , W.sub.L.sup.E) and the parameters of the decoder .sub.D=(W.sub.1.sup.D, . . . , W.sub.L.sup.D) define the parametrized encoder .sub.E(x)=.sub.L.sup.E(W.sub.L.sup.E . . . .sub.1.sup.E(W.sub.1.sup.Ex)) and the parametrized decoder .sub.D(x)=.sub.L.sup.D(W.sub.L.sup.D . . . .sub.1.sup.D(W.sub.1.sup.Dx)) with D layers, respectively. The latent space variables {u.sub.i}.sub.i=1.sup.M(310) represent a sort of compressed version of the set of training signals {x.sub.i}.sub.i=1.sup.M(20). The model training, i.e. estimation of {circumflex over ()}.sub.E and {circumflex over ()}.sub.D, can be performed similarly to the predictive sparse decomposition technique mentioned in the introduction, which leads to equation

[00003]$\begin{array}{cc}\left({\hat{}}_E,{\hat{}}_D\right)={\mathrm{argmin}}_{{}_E,{}_D}.Math.\underset{i=1}{\overset{M}{.Math.}}.Math.\frac{1}{2}.Math.{.Math.x_1-{}_{{}_D}\left(u_i\right).Math.}_2^2+{}_u.Math.{}_u\left(u_i\right)+\frac{}{2}.Math.{.Math.u_i-{}_{{}_E}\left(x_i\right).Math.}_2^2+{}_1.Math.{}_{{}_E}\left({}_E\right)+{}_2.Math.{}_{{}_D}\left({}_D\right),& \left(18\right)\end{array}$

where .sub..sub.E and .sub..sub.D(.) are technical constraints on the parameters of the encoder and the decoder, respectively. Assuming .fwdarw., this results into u.sub.i=.sub.E(x.sub.i) and leads to the generation of a broad family of auto-encoders in a form that is, in principle, known in the context of so-called deep learning techniques mentioned in the introduction and that can be described by equation

[00004]$\begin{array}{cc}\left({\hat{}}_E.Math.{\hat{}}_D\right)={\mathrm{argmin}}_{{}_E,{}_D}.Math.\underset{i=1}{\overset{M}{.Math.}}.Math.\frac{1}{2}.Math.{.Math.x_i-{}_{{}_D}\left({}_{{}_E}\left(x_i\right)\right).Math.}_2^2+{}_u.Math.{}_u\left({}_{{}_E}\left(x_i\right)\right)+{}_1.Math.{}_{{}_E}\left({}_E\right)+{}_2.Math.{}_{{}_D}\left({}_D\right).& \left(19\right)\end{array}$

[0068] In this context, one can use the ADMM approach mentioned in the introduction and known as such by a person skilled in the art, in combination with efficient optimizers such as ADAM, also mentioned in the introduction and known as such by a person skilled in the art, at each stage.

[0069] Alternatively, one can also consider even simpler approaches such as k-means applied to the subbands of the transform domain x to train the encoder-decoder pair. This alternative will be explored further below for reasons of demonstrating the power of the method according to the present invention, even if used in combination with such a simple model as k-means for determining the learnable priors, over the above described first and second cases wherein no priors, respectively only hand-crafted priors are used and wherein the decoder does not use the learnable priors and only the sampling operator is trained.

[0070] Once the encoder-decoder pair ({circumflex over ()}.sub.E, {circumflex over ()}.sub.D) is trained on the set of training signals {x.sub.i}.sub.i=1.sup.M, one can reformulate equation (4) using the parametrized regularization prior .sub.x(x) into cost functions described by the following equations

[00005]$\begin{array}{cc}{J}_1\left(u,x\right)=\frac{1}{2}.Math.{.Math.y-P_{}.Math..Math..Math.x.Math.}_2^2+\frac{}{2}.Math.{.Math.x-{}_{{}_E}\left(u\right).Math.}_2^2,& \left(20\right)\\ J_2\left(u,x\right)=\frac{1}{2}.Math.{.Math.u-{}_{{}_D}\left(x\right).Math.}_2^2+{}_u.Math.{}_u\left(u\right).& \left(21\right)\end{array}$

[0071] Applying the Nash equilibrium mentioned in the introduction and known by a person skilled in the art to equations (20) and (21) similarly to leads to a consensus of restrictions imposed by said cost functions hand J.sub.1 and J.sub.2. Said consensus of restrictions is defined as a fixed point (u*, x*) fulfilling equations

x*=argmin.sub.xJ.sub.1(u*,x*), (22)

u*=argmin.sub.uJ.sub.2(u,x*), (23)

this being similar in spirit to the above mentioned ADMM approach.

[0072] One can demonstrate that a solution to equations (22) and (23) is obtained, by skipping here the dependence of the estimate of the signal of interest on the set of sampling parameters , as an iterative decoder of a form described by equation

{circumflex over (x)}.sup.k=Fy(1FP.sub.){tilde over (x)}.sup.k, (24)

with an estimate of the signal of interest {tilde over (x)}.sup.k=({circumflex over (x)}.sup.k1, .sub.u) where (.) is a non-linear function depending on a set of sampling parameters .sub.u(.) and a transform domain F=(*P.sub..sup.TP.sub.).sup.+P.sub..sup.T, with .sup.+ denoting a pseudo-inverse, that can be further reduced for special forms of the transform domain .

[0073] Moreover, one can also consider, for equation (24), to use the LISTA implementation which is mentioned in the introduction as well as known by a person skilled in the art and which has a form of a deep network such as shown schematically in FIG. 4. The latter schematically illustrates an implementation of an iterative decoder in a form of a multi-layer network with fixed blocks of model encoding and decoding as well as common reconstruction blocks A and B chosen to be fixed for illustration purposes. For the common reconstruction blocks A and B, it is assumed that A=F and B=1FP.sub.. In a more general case, the model encoder-decoder blocks as well as blocks B can be updated at each layer according to the minimization of the overall cost function that will increase the quality of reconstruction and minimize the number of layers needed in comparison with the above mentioned example of a fixed block network decoder.

K-Means Based Model Encoder-Decoder with Subbands Partition

[0074] In this section, it will be presented how the model training may be implemented in the context of the proposed joint learning of the sampling operator and the decoder by use of a k-means model applied to the subbands of the transform domain x to train the encoder-decoder pair. The procedure again comprises two stages, namely training/learning and testing/sampling.

[0075] The training/learning stage may be described by the following procedure:

[0076] 1. One assumes some complex transform domain consisting of =(.sub.re,.sub.im) that leads to a representation described by equation

x=.sub.rex+j.sub.imx=x.sub.re+jx.sub.im. (25)

such as Fourier transform.

[0077] 2. The Fourier spectrum is split into L subbands. The splitting might be overlapping or non-overlapping. The splitting is done such as to ensure that the magnitude of the Fourier transform domain components of each subband

[00006]$e\left(l\right)=\frac{1}{V_l}.Math.\sqrt{\underset{b=1}{\overset{V_l}{.Math.}}.Math.{.Math..Math..Math.x_b\left(l\right).Math.}^2},l\left(1,\mathrm{.Math.}.Math.,L\right)$

is the same, where V.sub.l is the number of samples in each subband l.

[0078] 3. Given a set of training signals X=(x.sub.1, . . . , x.sub.M), the latter is split into L subbands as (X)=((X.sub.1), (X.sub.L)).

[0079] 4. Thereafter, two codebooks for the real C.sub.re and imaginary C.sub.im parts are generated using the k-means algorithm which is well known to a person skilled in the art. One can generate the codebooks independently or jointly.

[0080] 5. The sampling operator training following thereafter consists of two parts that lead to minimization of the overall reconstruction error. The first part aims at minimizing the reconstruction error from the direct observation signal y in the set of sampling parameters . It naturally leads to selection of coefficients with the largest magnitudes in the Fourier spectrum X of the set of training signals. The second part seeks minimizing the error between the set of centroids and the samples in the set of sampling parameters . At the same time, it is important to note that the set of sampling parameters should be chosen in such a way to ensure the maximum distinguishability of codewords in the codebooks. Otherwise, if the distances to several codewords in the sampling points of the set of sampling parameters are very small, the codewords are not distinguishable and any random perturbation might flip them. The result of this procedure will lead to a set of sampling parameters ensuring the most informative and accurate selection of the codewords in the trained codebooks C.sub.re and C.sub.im. The considered trained representation also closely resembles the compressed representation of data in the trained codebooks.

[0081] The testing/sampling stage may be described by the following procedure:

[0082] Given an observation signal which may be represented as a vector y=P.sub.x=(P.sub.1x,P.sub.2x, . . . , P.sub.Lx) as a subband based sampling, the decoder first finds the closest representatives in the codebooks and C.sub. for each subband l by use of equation

c.sub.re.sub.l()=argmin.sub.1jK.sub.ly.sub.re.sub.lP.sub..sub.lc.sub.re.sub.l(j).sup.2.sub.2

c.sub.im.sub.l()=argmin.sub.1jK.sub.ly.sub.im.sub.lP.sub..sub.lc.sub.im.sub.l(j).sup.2.sub.2, (26)

where K.sub.l denotes the number of codewords in each subband l(1, . . . , L).

[0083] One can also include the magnitude and phase components to improve the overall system performance, such as described by equation

[00007]$\begin{array}{cc}\left(c_{{\mathrm{re}}_l}\left(\hat{J}\right),c_{{\mathrm{im}}_l}\left(\hat{J}\right)\right)={\mathrm{argmin}}_{1j{K}_l}.Math.{.Math.y_{{\mathrm{re}}_l}-P_{{}_l}.Math.c_{{\mathrm{re}}_l}\left(j\right).Math.}_2^2+{.Math.y_{{\mathrm{im}}_l}-P_{{}_l}.Math.c_{{\mathrm{im}}_l}\left(j\right).Math.}_2^2+{}_1.Math.{.Math..Math.y_l.Math.-.Math.P_{}\left(c_{{\mathrm{re}}_l}\left(j\right)+{\mathrm{jc}}_{{\mathrm{im}}_l}\left(j\right)\right).Math..Math.}_2^2+{}_2.Math.{.Math.\arg \left(y_l\right)-\arg \left(P_{}\left(c_{{\mathrm{re}}_l}\left(j\right)+{\mathrm{jc}}_{{\mathrm{im}}_l}\left(j\right)\right)\right).Math.}_{2.Math.}& \left(27\right)\end{array}$

where |.| and arg(.) denote the magnitude and phase of a complex vector, respectively, ..sub.2 is a mean angular difference normalized in the range (0,2), and .sub.1 and .sub.2are Lagrangian multipliers.

[0084] The reconstruction is based on the decoder producing an estimate of the signal of interest according to equation

[00008]$\begin{array}{cc}\hat{x}={}^{*}.Math.\left\{\begin{array}{c}\left(P_{{}_1}^T.Math.y^1+P_{{}_1^C}^T.Math.P_{{}_1^C}\left(c_{{\mathrm{re}}_1}\left(j\right)+{\mathrm{jc}}_{{\mathrm{im}}_1}\left(j\right)\right)\right),\\ \mathrm{.Math.},\\ \left(P_{{}_L}^T.Math.y^L+P_{{}_L^C}^T.Math.P_{{}_L^C}\left(c_{{\mathrm{re}}_L}\left(j\right)+{\mathrm{jc}}_{{\mathrm{im}}_L}\left(j\right)\right)\right),\end{array}\right\}& \left(28\right)\end{array}$

where {.} denotes the concatenation of vectors. The main difference of this decoder as compared to the linear decoder disclosed in the article Learning-based subsampling and the document US 2017/0109650 consists in the presence, in the decoder according to the present invention, of injected learnable priors coming from the learned K-means codebooks. Another difference between the decoders is the presence, in the decoder according to the present invention, of multiple subbands optimized for optimal sampling.

Choice of the Cost Function

[0085] In general, the cost function used in the sampling method according to the present invention at the cost minimization step (180) for determining the set of sampling and decoding parameters .sub.g during the joint training of the sampling operator P.sub.(120) and of the decoder g.sub.g(.) (140) is chosen depending on the targeted application. In particular, the cost function may comprise a metric of closeness between two signals of interest x(60) in the direct domain or latent space for recognition of the signal of interest x(60) or may comprise a metric of closeness between the targeted and estimated labels for classification, recognition, authentication, identification or forensic analysis of the signal of interest x(60). Furthermore, the cost function used in the cost minimization step (180) for determining the set of sampling and decoding parameters .sub.g during the joint training of the sampling operator p.sub.(120) and of the decoder g.sub.g(.)(140) may optionally contain an additional constraint on the latent space ensuring rate-efficient representation of model parameters for a joint minimization of the sampling rate in the sampling operator P.sub. and of the compression rate in the model latent space for a desired reconstruction distortion used for joint sampling-compression applications. These features are particularly interesting in applications of the present method related to security, like identification or authentication.

Modified Cost Function for Classification and/or Recognition

[0086] For applying the method according to the present invention to the above mentioned classification problems, in particular to classification and/or recognition of the signal of interest x based on the observation signal y, the training stage is formulated as an optimization problem according to equation

({circumflex over ()},{circumflex over ()}.sub.g)=argmin.sub.,.sub.g.sub.p=1.sup.P.sub.i=1.sup.M.sub.m=1.sup.CL(l(m),g.sub.g(P.sub.(x.sub.i+z.sub.p)))+.sub.1.sub.()+.sub.2.sub.g(.sub.g), (29)

wherein the decoder g.sub.g(.) is a classifier, i.e. g.sub.g:.sup.n.fwdarw.{1, . . . , C} and l(m) is an encoded class label for the class g.sub.gL.sup.n.fwdarw.{1, . . . , C}. L(.) denotes a corresponding cost function typically used 1.sup.1:

[0087] in classification problems. For reasons of clarity, it is added here that this description refers to the classification problem as a general problem. The classification problem can be split into two problems, namely if there are only two hypothesisit is referred to as authentication problem, and if there are several hypothesisit is referred to as recognition or identification problem. For the sake of completeness, it shall be noted that the terms recognition and identification are sometimes used with different signification in research communities. In the present description, identification shall signify a classification problem having M+1 hypothesis, i. e. including a reject option in case the probe does not fit to any of the known classes, whereas recognition shall signify a classification problem having M hypothesis and the probe should be classified as one of the known classes without a reject option.

Modified Cost Function for Authentication

[0088] In order to allow for authentication by the method according to the present invention, the training stage is formulated as an optimization problem according to equation

({circumflex over ()},{circumflex over ()}.sub.g)=argmin.sub.,.sub.g.sub.p=1.sup.P.sub.i=1.sup.M.sub.m=0.sup.1L(l(m),g.sub.g(P.sub.(x.sub.i+z.sub.p)))+.sub.1.sub.()+.sub.2.sub.g(.sub.g), (30)

wherein the decoder g.sub.g(.) is a binary authentication classifier, i.e. g.sub.g:.sup.n.fwdarw.{0,1} and l(m) is an encoded class label. L(.) a gain denotes a corresponding cost function.

[0089] Both approaches for classification based on the recognition and on the authentication formulations are shown in FIG. 6 which schematically illustrates a joint training of the sampling operator and the decoder for a generic classification problem. The cost function is optimized at the training stage to ensure the desired output of the decoder under the constraint on the sampling operator. In contrast to FIG. 2, where the output of the decoder is an estimate of signal {circumflex over (x)}.sub., the decoder in the classification setup outputs an estimate of class label {circumflex over (l)}.sub.(m) for the specified class m={1, . . . , C}. The class-dependent set of training signals {x.sub.i(m)}.sub.i=1.sup.M, m={1, . . . , C} is labelled accordingly in FIG. 6.

[0090] FIG. 7 shows a first embodiment of a multi-stage sampling-decoding framework wherein the signal of interest x(60) is first sampled by the sampling operator P.sub..sub.1(600) producing a sampled signal y.sub.1(610) that might represent a low-resolution sampled version of the signal of interest x. The first decoder (620) produces an output (630) that might correspond to the detection region of some desired features that in turn might be sampled with different rates. One might consider a face detector from a low-resolution image as a specific example of this general approach. Finally, the sampling operator P.sub..sub.2(640) then samples the trained regions to produce the observation signal y.sub.2 whereas the second decoder (660) performs the recognition/detection by producing an estimate (520). This example is mentioned mainly for illustrative purposes, given that one can apply different optimization strategies to this example in various applications which cannot all be mentioned here, the optimization moreover depending on the given energy consumption, exposure time, memory, complexity, security or privacy considerations for a given application. One can also envision application of the above described multi-stage sampling-decoding scenario to the case where different sensors perform sampling of different regions. Additionally, this scenario is illustrated in FIG. 7 for two stages only, whereas in practice multi-stage systems having several stages might be deployed to fulfill the targeted requirements. The overall cost functions at each stage might have different objectives, too.

[0091] FIG. 8 shows a second embodiment of a multi-stage sampling-decoding framework which extends the example of FIG. 7 to a transform domain (670) with several sampling operators, wherein first the sampling using the sampling operator P.sub..sub.1(600) is performed in some domain (680) and the sampling using the sampling operator P.sub..sub.2(640) in another domain (630), based on the result of the first stage decoder 1 (620). For example, the first transform domain (680) may be chosen to be the identity operator =1, such that the first sampling using the sampling operator P.sub..sub.1(600) is performed in the spatial domain instead of in some transform spectral domain, whereas the next sampling using the sampling operator P.sub..sub.2(640) indeed is performed in a spectral domain (630).

[0092] The proposed method can be used in several kind of imaging applications. Without being exhaustive, only a few of a high number of potential applications will be considered explicitly in the following by way of example, bearing in mind that analoguous setups may easily be transposed by a person skilled in the art to applications not explictly mentioned in the following.

[0093] Application of the method according to the present invention to medical imaging may include, but is not limited to, magnetic resonance imaging (MRI), computerized tomography (CT) and ultrasound imaging. In all these applications, the image is sampled in some transform domain represented by physics of the corresponding imaging problem.

[0094] Non-destructive testing and imaging include applications where the sampling is performed in some transform domain Wand one is interested to visualize the scene under investigation or to make some decision on that.

[0095] Sparse sensor arrays include numerous applications in radar, radioastronomy, acoustic arrays, surveillance systems and remote sensing applications.

[0096] Applications related to CCD and CMOS sensors include all professional, multimedia and specialized applications where use of huge amounts of pixels is limited or not desirable due to some technical limitations that include but are not limited by energy consumption, memory, communication burden, etc. Additionally, it might be interesting to use systems with larger size of sensors to increase the sensitivity without a big sacrifice in the resolution to cope with the photon noise. Moreover, one can also develop efficient sampling schemes adapted to the compression requirements of Big Data applications. Finally, one can also significantly benefit from the proposed sampling and the corresponding demosaicking as a kind of decoder in color imaging applications.

[0097] The invention is also of interest for classification and recognition applications when some object or human should be recognized from a limited number of samples. This situation is of interest not only from technical perspectives but also from security and privacy considerations. In particular, it is possible to construct a working system if the number of samples is just sufficient for performing a reliable identification or authentication, but is insufficient to reliably reconstruct the whole object or biometric parameters. Such a sampling is also of interest in order to develop new countermeasures against the adversarial attacks in machine learning. Additionally, it is possible to randomize the sets of sampling parameters .(fusing secret keys jointly with sets of training signals. Additionally, it is possible to introduce permutation of acquired samples by way of an operator Q:.sup.n.fwdarw..sup.n in equation (1), which insures key-based permutations. In this way, one can create a non-differentiable operator preventing adversarial learning and attacks.

[0098] The proposed joint training of the sampling operator and of the decoder is also of interest for anticounterfeiting applications, in which case the sampling operator and the decoder are trained on a set of training signals consisting of authentic objects/biometrics and of fakes with the above mentioned binary authentication cost function. In this manner, the found sampling regions allow to quickly and reliably establish the authenticity of objects or identity of humans without need of providing the whole range sampling. Additionally, this leads to more efficient implementations and increased reliability.

[0099] Finally, the proposed invention has a wide range of applications for both image reconstruction and classification problems in case the signals of interest represent high-dimensional spatio-temporal tensors. In these cases, the proposed sampling method might be especially beneficial in comparison to prior art sampling methods. FIG. 9 illustrates application of multi-stage sampling to spatio-temporal signals (800). The signal of interest represented by tensor (800) is sampled in a first step in the spatio-temporal group x.sup.(1)(810) by the sampling operator P.sub..sub.1(820), producing observation signal y.sub.1(830). The decoder 1 (840) produces an estimate for the sampling operator P.sub..sub.2(870) applied in the next step of sampling to the spatio-temporal group (x.sup.(2), . . . , x.sup.(J))(860), resulting in y.sub.2(880) and the decoder 2 (890) produces an estimate of class label (900). In general, the resulting decoder can also produce an estimate of the signal of interest {circumflex over (x)} depending on the targeted cost function.

[0100] In order to illustrate the results that may be achieved by application of the method according to the present invention in comparison to prior art methods, FIG. 10 presents the results of sampling and reconstruction for one medical (MRI) and two astronomical datasets of images denoted as D1, D2 and D3, respectively, for the Euclidean cost function for different sampling rates. All images have been cropped to a size of 256256 pixels. The results of two embodiments of the method according to the present invention, both based on learnable priors, are compared to the prior art on-fly S-best sampling method disclosed in the article Robust image restoration matched with adaptive aperture formation in radar imaging systems with sparse antenna arrays by I. Prudyus et al., denoted in FIG. 10 as S-best, and to the prior art sampling method disclosed in the article Learning-based subsampling as well as in the document US 2017/0109650 by V. Cevher et al., denoted in FIG. 10 as Average. The first embodiment of the method according to the present invention, denoted in FIG. 10 as Proposed K-means prior, is based on the above described K-means learning of codebooks adapted to the Fourier spectral subbands of test image datasets having the same number of centroids per subband that can be itself a factor of additional optimization, leading to an additional performance enhancement. The second embodiment of the method according to the present invention, denoted in FIG. 10 as Proposed Deep image prior, is based on the above mentioned decoder learning which, in the embodiment the results of which are shown in FIG. 10, uses deep learning, the latter being mentioned in the introduction and being as such known in prior art. The second embodiment is supposed to illustrate the potential use of deep learning methods within the limited training on external sets of training signals. In all cases, the traditional mean square error (MSE) computed between the original and reconstructed images was used to estimate the quality of reconstruction. Additionally, to reflect the aspect of human subjective perception of reconstructed image quality, the well-known SSIM metric known by a person skilled in the art was used. As is clearly visible in FIG. 10, both embodiments of the method according to the present invention comprising joint training of the sampling operator and decoder, even if using a very simple k-means model of trainable priors, significantly outperforms all tested methods for both the MSE and SSIM metrics. It should be noted here that the particular performance enhancement investigated on broader sets of images ranging from real life photos to specialized domain imagery depends on the uniformity of the data within the training datasets and its correspondence to the signals of interest. This highlights that proper training and use of learnable priors is beneficial and may greatly improve image sampling methods.

[0101] For further illustration, FIG. 11 illustrates some examples of images that were tested, i.e. sampled, in order to obtain the results shown in FIG. 10. As suggested by the objective MSE and SSIM measures of FIG. 10, the quality of the images illustrated in FIG. 11, sampled and reconstructed for very small sampling rates, is significantly better for both embodiments of the method according to the present invention in comparison to the prior art S-best and Average techniques. The quality increases with the increase of the sampling rate, where all methods produce quite similar results.

[0102] Finally, it shall be noted here that the present invention also pertains to computer program means stored in a computer readable medium adapted to implement the above described method as well as to a device equipped with such computer program means. The device preferably is chosen from the group comprising a mobile phone, in particular a smart phone equipped with a camera, a digital photo apparatus, a digital video camera, a scanning device, a tablet or personal computer, a server.

[0103] In light of the above description of the method of signal sampling according to the present invention, its advantages are clear. In particular, due to the joint adaptation of the sampling operator and decoder to the training data, the method of signal sampling according to the present invention has a higher overall performance under a targeted cost function for a signal of interest acquired with a finite number of samples. This allows to reduce the number of samples in many applications that have constraints in terms of technology, cost, memory and exposure/observation time as well as in terms of health safety, security or privacy, etc. by simultaneously ensuring high quality of decoding. Furthermore, the joint training of the sampling operator and of the decoder allows, in contrast to prior art sampling methods, to avoid use of hand-crafted decoders, to provide improved decoders most adapted to the corresponding sampling operators, as well as to exploit training signals to its full power. Moreover, the proposed method allows to make a decision about a class of the signal of interest in recognition and identification applications and the training procedure of the sampling operator may be improved by taking into account in more precise manner as compared to prior art the predetermined properties of specific types of signals, for example of the kind of images, to be sampled and reconstructed, respectively classified.