EP4264553A1

EP4264553A1 - Method for reconstructing an image based on compressed acquisition of measurements

Info

Publication number: EP4264553A1
Application number: EP21823852.5A
Authority: EP
Inventors: Antoine Souloumiac
Original assignee: Commissariat a lEnergie Atomique CEA; Commissariat a lEnergie Atomique et aux Energies Alternatives CEA
Current assignee: Commissariat a lEnergie Atomique et aux Energies Alternatives CEA
Priority date: 2020-12-17
Filing date: 2021-12-01
Publication date: 2023-10-25
Also published as: WO2022128465A1; FR3118251A1; FR3118251B1

Abstract

Method for reconstructing an image I of an object comprising the steps of: - acquiring (201), in a vector m, a number p of object imaging measurements; - determining (202) a first matrix A^T equal to the product of a projection matrix P defining a set of projection vectors and of a dictionary matrix D formed of the components of an image decomposition basis, - seeking (203) a parsimonious solution x, of size N strictly greater than d, to the underdetermined linear system A^Tx=m - reconstructing (204) the image I through the product of the dictionary matrix D and of the vector x, which is the solution to said system.

Description

DESCRIPTION

Title of the invention: Method for reconstructing an image from a compressed acquisition of measurements

The invention relates to the field of imaging, for example imaging by X-rays or by magnetic resonance or any other imaging method making it possible to reconstruct a 2D or 3D image of an object. The invention applies in particular to the field of medical imaging or astronomical imaging.

The invention also relates to the field of compressed acquisition or "compressed sensing" in English which aims to reconstruct a signal from a limited number of acquisitions of measurements of this signal.

[0003] In order to be able to reconstruct an image of a given object with a precise resolution, this requires carrying out a large number of acquisitions by means of an imaging device, according to different viewing angles.

[0004] For example, in the case of X-ray imagers consisting of an X-ray source and a set of detectors, the source and the detectors are positioned on either side of the object to be imaged. and perform a measurement of the X-ray absorptions of a voxel of the object along a direction defined between the source and the detectors. To reconstruct an image of the object, the source and the detectors rotate around the object in order to modify, for each acquisition, the direction, that is to say the angle of view.

[0005] The greater the number of acquisitions from different viewing angles, the better the resolution of the reconstructed image. However, each acquisition has a cost in terms of acquisition time by the imaging device, energy consumed by the device or even X-ray dose absorbed by the patient (which must be limited to limit the side effects harmful effects of ionizing radiation).

[0006] There is therefore a need to reduce the number of acquisitions while maintaining satisfactory image quality. [0007] The known techniques of compressed acquisition ("compressed sensing") are based on the search for a parsimonious solution to an under-determined linear system.

[0008] An algorithmic difficulty resides in the calculation of this parsimonious solution. Various state-of-the-art solutions have been proposed to solve this problem: "matching pursuit" or "orthogonal matching pursuit" type algorithms, thresholding techniques or even involving the L ₁ norm. These solutions are by example described in references [1], [2], [3] and [4],

[0009] All of these methods have drawbacks related to problems of convergence, accuracy of the solution or complexity of implementation.

The invention proposes a new compressed acquisition method based on the search for linear combinations of maximum kurtosis to determine a parsimonious solution to the aforementioned underdetermined linear system.

[001 1] The compressed acquisition method according to the invention makes it possible to reduce the number of acquisitions of an imaging device in order to reconstruct an image without loss of quality.

The subject of the invention is a method for reconstructing an image I of an object, the image I having a number d of pixels, the method comprising the steps of:

- Acquire, in a vector m, a number p, strictly less than d, of imaging measurements of the object, by means of an imaging device,

- Determine a first matrix ^AT equal to the product of a projection matrix P defining a set of projection vectors and a dictionary matrix D made up of the components of an image decomposition base,

- Find a sparse solution x, of dimension N strictly greater than d, to the underdetermined linear system A ^T x=m in: i. Determining a solution x _M c of minimum quadratic norm of said linear system and normalizing the vector x _M c ii. Determining a second matrix whose rows form an orthonormal basis of the complementary subspace orthogonal to the rows of the first matrix A ^T , iii. Finding a linear combination x of the vector x _M c and the rows of the second matrix having maximum kurtosis.

- Reconstruct the image I by the product of the dictionary matrix D and the vector x, solution of the said system.

According to a particular aspect of the invention, said linear combination is obtained by searching, for at least one row of the second matrix , a linear combination of a pair of components formed by the vector x _M c and this line, having a maximum kurtosis.

[0014] According to a particular aspect of the invention, said linear combination is obtained by searching for a Givens rotation matrix, defined by a rotation angle θ , for which the product of said pair of components and of this Givens rotation matrix exhibits a maximum kurtosis for one of its two components or exhibits a maximum difference between the kurtosis of the two components.

According to a particular aspect of the invention, the search for a linear combination of maximum kurtosis is iterated for several rows of the second matrix

According to a particular aspect of the invention, the dictionary matrix D consists of a base of elementary signals, for example a base of wavelet functions, or is generated by automatic learning from learning images.

According to a particular aspect of the invention, the vector m is acquired by means of an X-ray detector and being made up of the concatenation of several measurement vectors, each measurement vector comprising several measurements carried out by the detector for a different angle of view of the object.

According to a particular aspect of the invention, the projection matrix P contains projection vectors respectively defining directions associated with an angle of view.

According to a particular aspect of the invention, the image I is reconstructed in two dimensions or in three dimensions. The invention also relates to a system for reconstructing an image I of an object, the image I having a number d of pixels, the system comprising an imaging device for acquiring, in a vector m , a number p, strictly less than d, of object imaging measurements and a processing unit configured for:

- Find a sparse solution x, of dimension N strictly greater than d, to the underdetermined linear system A ^T x=m in: i. Determining a solution x _M c of minimum quadratic norm of said linear system and normalizing the vector x _M c ii. Determining a second matrix whose rows form an orthonormal basis of the complementary subspace orthogonal to the rows of the first matrix ^AT , iii. Finding a linear combination x of the vector x _M c and the rows of the second matrix A^ having maximum kurtosis.

According to a particular aspect of the invention, the processing unit is configured to execute the steps of the image reconstruction method according to the invention.

According to a particular aspect of the invention, the imaging device is an X-ray detector or a magnetic resonance imaging unit.

According to a particular aspect of the invention, the system comprises a display device for displaying the reconstructed image I.

Other characteristics and advantages of the present invention will appear better on reading the following description in relation to the following appended drawings. [0025] Figure 1 shows a diagram illustrating the principle of an X-ray imager,

[0026] FIG. 2 represents a flowchart detailing the steps for implementing an image reconstruction method according to the invention,

[0027] FIG. 3 represents a flowchart detailing the steps for implementing the search for a sparse solution of an underdetermined linear system according to an embodiment of the invention.

[0028] Figure 1 schematically illustrates the operating principle of an X-ray imager or X-ray scanner comprising an SRX X-ray source, for example an X-ray tube and a DRX detection device consisting of several elementary detectors for example arranged linearly or in the form of a matrix. The object to be imaged is placed between the SRX source and the DRX detection device and is irradiated by RX X-rays emitted by the source towards the detectors in a direction predetermined by an angle θ. Each pixel of the object is characterized by an X-ray absorption value and it is sought to reconstruct the image, denoted l(x,y), of these values. The measurement made by a detector corresponds to the combination of these absorption values along the ray connecting the SRX source to the detector.

The elementary detectors of the DRX device therefore measure the projection of an image l(x, y) on these detectors in the direction defined by the angle Q.

[0030] To completely image the object, the assembly consisting of the SRX source and the DRX detectors pivots around the object in order to carry out several successive acquisitions at different angles Q.

In the field of medical imaging, the object to be imaged corresponds to a thin anatomical slice or section. The SRX source and the DRX detectors are constrained to be positioned in the plane of the anatomical slice and rotate around the slice.

The two-dimensional image of the anatomical slice is reconstructed from the different projections measured by the detectors for different angles Q.

The same principle of acquisition can be achieved in 3D tomography by generalizing this principle to a body volume. The SRX source and detectors XRDs pivot in this case around the volume to be imaged according to optimized trajectories ensuring the most complete coverage possible of the different angles of view. The directions of projection are, in the case of 3D imaging, defined by two angles which determine the point of view of the imaging device with respect to the volume to be imaged.

One objective of the invention is to provide a method making it possible to reconstruct the image I (in 2D or in 3D) from a reduced number of acquisitions.

Figure 2 details the steps for implementing the image reconstruction method according to one embodiment of the invention.

The first step 201 consists in acquiring several imaging measurements by means of an imaging device.

[0037] In the case where the imaging device used is an X-ray scanner of the type described in FIG. 1, each measurement consists of a vector of N _d values corresponding to the number of elementary detectors of the ray detection device X XRD. N _a measurements are carried out for various different viewing angles θ as described above.

The acquired measurements are digitized and concatenated in a vector m of size p=N _d ×N _a . One objective of the invention is to reconstruct an image I of the object to be imaged, the image having a number d of pixels strictly greater than the dimension p of the vector m.

Without departing from the scope of the invention, any other imaging device capable of producing a vector of imaging measurements m of an object, according to different projections, can be used as a replacement for the X-ray scanner. In particular, a magnetic resonance imaging (MRI) unit can be used to perform the measurement step 201. In this case, the projection angle variation described above for the X-ray scanner is replaced by a variation of a magnetic field created by currents passing through coils called field gradient coils. This variable field is added to the permanent field, fixed in time and useful space, created by a superconducting magnet.

These field variations are modeled as a trajectory in an abstract space called "k-space" in English. An MRI sequence consists of a certain trajectory to sample this space, the more samples the better the image. The time required to acquire an image is highly variable: from the order of one second for functional MRI to several tens of minutes for a T2 image, but in all cases it is useful to shorten this time without degrading the quality. image.

[0041] For each given position in “k-space”, the MRI signal is collected by radio frequency antennas placed close to the patient inside the superconducting magnet.

[0042] Thus, for a magnetic resonance imaging application, the vector m corresponds to different measurements taken by the antennas of the MRI unit, these measurements being associated with a given position in "k-space" space. .

Returning to the initial example of an X-ray scanner, if we note P the projection matrix which contains the vectors associated with the directions given by the different viewing angles 9, we can write:

[0044] m= P.l,

Indeed, the vector m corresponds to the result of the projection of the desired image I on the projection base defined by the matrix P.

Furthermore, the image I can be expressed as a sparse linear combination of reference images taken from a dictionary D defining an image decomposition base:

[0047] l= D.x with D a matrix of dimensions d by N (N>d) and x a sparse vector of size N, that is to say comprising a large proportion of zero values.

The matrix D contains the components of a signal decomposition base, for example a decomposition base into wavelets or into ridgelet functions or into curvelet functions or any other signal base which makes it possible to decompose the image I .

In the case where the image I to be acquired belongs to a particular application domain, the dictionary matrix D can be obtained at the end of automatic dictionary learning from learning images, by example via the method described in reference [5]. [0050] From the previous relations, we can therefore write m=PDx=A ^T .x.

In a step 202, the matrix A ^T = PD (A ^T ap rows and N columns) is determined from the projection matrix P and the dictionary D, these two matrices P and D being determined as a function of the intended application.

In a step 203, the most parsimonious possible solution x to the underdetermined linear system A ^T .x =m is then sought.

Finally, the image l=D.x is reconstructed (step 204) from the solution x determined in step 203.

FIG. 3 schematizes the detail of step 203 making it possible to determine the solution x.

The calculation of a solution x as parsimonious as possible of the preceding underdetermined linear system can be reduced to the maximization of a kurtosis in the following way. We can always consider that the parsimonious solution x is equal to the sum of:

- a solution of minimum quadratic norm x _MC = AA ^T A')~ ¹ m (this solution is a linear combination of the rows of the matrix A ^T and is, in general, not parsimonious) and

- of a vector orthogonal to the rows of the matrix A ^T that we note x _± . If we note A a matrix with Np rows, and N columns and whose rows form an orthonormal basis of the complementary subspace orthogonal to the p rows of A ^T , then is a linear combination of the columns of 4 _± .

We can verify that x = x _MC + x _± is indeed a solution of the system A ^T x = m whatever x _± since

Thus, the search for the most parsimonious possible solution x is equivalent to the search for a linear combination of x _MC and the rows of A which is the most parsimonious possible.

This amounts to finding the linear combination of the rows of maximum kurtosis of the matrix X equal to [0059] Indeed, a very high kurtosis characterizes a random variable very close to 0 with a high probability and non-zero with a very low probability, in other words a parsimonious random variable.

In other words, step 203 of the method according to the invention is carried out by means of the sub-steps described in FIG. 3:

- Step 301: determination of the solution of minimum quadratic norm XMC ⁼ A(A ^T A') ^-1 m

- Step 302, determination of the matrix A^ whose rows form an orthonormal basis of the complementary subspace orthogonal to the rows of the matrix A ^T .

- Step 303: search for a linear combination of the vector x _MC and the rows of the matrix which exhibits maximum kurtosis.

The step 302 does not depend on the measurement vector m and can be performed initially once and for all and not for each image to be reconstructed.

The vector x _MC is normalized so as to obtain a matrix X with orthonormal rows.

We define the empirical kurtosis of a row (x ₁ ,x ₂ , ...,x _n , ...x _N ) of the matrix X by the formula This formula is simplified because of the rows orthonormal because we have

[0064] The simplified formula for empirical kurtosis becomes

The search for a maximum value of the empirical kurtosis K _x is equivalent to the search for a maximum value of the simplified quantity

Step 303 is performed by determining different linear combinations of the first row of the matrix X with the other rows.

In other words, step 303 is performed by determining a square rotation matrix U, of dimensions (N-p+1)x(N-p+1) and by multiplying X on the left by this matrix U, or U x X. The matrix X has dimension N-p+1 by N. We notice that it is advantageous to multiply on the left only by orthogonal matrices (matrices U which satisfy U ^T U = UU ^T = I _N-P+1 = the identity matrix of dimension (/V - p + 1) x (/ V - p + 1)) because they preserve the normalization of the lines that is to say the validity of the simplified formula of the empirical kurtosis.

The orthogonal matrix U such that the first row of U x X has maximum kurtosis cannot be calculated in a single step; it is constructed iteratively as a product (the product matrices of orthogonal matrices are orthogonal) of particular elementary orthogonal matrices called Givens rotations (which correspond to plane rotations) denoted G(i,j,0) and j ≤ N — p + 1 and 0 ≤ θ ≤ π.

Note that the product G(i,j, 0) x X = X' replaces the rows X(i,: ) and X(j,: ) of X by the rows and X'(j, : ) equal to the following linear combinations:

[0070] X'(i,: ) = cos θ X(i,: ) — sin θ X(j,: ) and X'(j,: ) = sin θ X(i,: ) + cos θ X( j, : ) .

[0071]X(i,:) denotes the line of index i of the matrix X. The term “:” is used to designate all the values of the index.

The step 303 of the method thus consists in optimizing the angle of rotation 0 for several pairs (l,j) of rows of the matrix X in order to maximize the kurtosis of the first row.

The use of Givens rotations thus decomposes the global optimization problem into a series of 2-dimensional optimization problems that are easier to perform.

Two variant embodiments can be envisaged for determining the angle of rotation 0 which makes it possible to obtain a maximum empirical kurtosis for the first row of the matrix.

[0075] A first embodiment variant consists in finding the angle of rotation

0 which maximizes the difference in kurtosis between the two lines of the couple. Indeed, in maximizing this deviation, this makes it possible to obtain an approximately maximum kurtosis for the highest kurtosis line.

Let X(1, : ) = (x _1, ...x _n , ... x _N ) and X(j, : ) = (y _1, ...y _n , ...y _N ') two rows of the matrix X for which we seek a linear combination of maximum kurtosis.

We define the two linear combinations u and v by:

[0079] Assuming the rows X(1, : ) and X(j, : ) orthonormal, we obtain orthonormal vectors u and v, their respective empirical kurtosis are defined by

The angle of rotation 0 which maximizes the difference of the empirical kurtosis K _u -K _v is then sought. In the following it is considered by convention that “u” corresponds to the line of highest kurtosis and “v” corresponds to the line of lowest kurtosis. Thus, we seek to maximize the difference of the empirical kurtosis K _u - K _v .

From expression (1), the values of cos 2θ and sin 2θ are determined which maximize this expression, then the values of cos θ and sin θ are deduced therefrom and finally the value of the angle of rotation θ which allows to obtain a maximum empirical kurtosis for the row u resulting from the rotation of the rows X(1, : ) and X(j, : ).

An advantage of this first variant is its low numerical resolution complexity. A second embodiment variant consists in directly seeking the angle of rotation θ which makes it possible to maximize the empirical kurtosis K _u of the line u which is equal to:

[0087]

The expression (2) corresponds to a linear and quadratic criterion under quadratic constraint. The angle 0 which maximizes this expression is obtained by numerical resolution.

One of the two variants described above is therefore applied to each pair of rows (1,j) of the matrix X according to a predetermined sequence. Different line scanning structures can be envisaged for this purpose.

For example, an iteration loop is performed by varying j from 2 to p+1 until convergence.

At each iteration (numbered by the index “k”), the products of the Givens rotation matrices obtained are accumulated in an orthogonal matrix U _k .

The iterative procedure is organized as follows.

U ₀ =I _N-P+1 and X ₀ =X are initialized, with I _N-P+1 the identity matrix of dimension (N - p + 1) x (N - p + 1).

Then at each iteration, for each determined rotation matrix, we calculate

The orthogonality of the Givens rotations (that is to say the property G (i,j, θ) ^T x G(i,j, θ) = I _N-P+1 ) guarantees the stability of the produces U _k x X _k during the iterations, indeed

We therefore have X _k =U _k x X throughout the optimization procedure. Once convergence is reached after a certain number K of iterations, the first line of is as sparse as possible, the vector x, sparse and system fix

The invention can be implemented by means of an imaging device of the type described in FIG. 1 coupled to a processing unit capable of executing the steps of the image reconstruction method, according to the invention, from the measurements m made by the imaging device.

Thus, the invention can be implemented as a computer program comprising instructions for its execution. The computer program may be recorded on a processor-readable recording medium.

[0099] Reference to a computer program which, when executed, performs any of the functions previously described, is not limited to an application program running on a single host computer. Rather, the terms computer program and software are used herein in a general sense to refer to any type of computer code (e.g., application software, firmware, microcode, or other form of computer instruction) that can be used to program one or more processors to implement aspects of the techniques described herein. The IT means or resources can in particular be distributed (“Cloud computing”), possibly using peer-to-peer technologies. The software code may be executed on any suitable processor (e.g., a microprocessor) or processor core or set of processors, whether provided in a single computing device or distributed among multiple computing devices (e.g. example as possibly accessible in the environment of the device). The executable code of each program allowing the programmable device to implement the processes according to the invention can be stored, for example, in the hard disk or in ROM. In general, the program or programs can be loaded into one of the storage means of the device before being executed. The central unit can control and direct the execution of the instructions or portions of software code of the program or programs according to the invention, instructions which are stored in the hard disk or in the ROM or else in the other aforementioned storage elements. [01001 References

[0101] [1] D.L. Donoho, "Compressed sensing," in IEEE Transactions on Information Theory, vol. 52, no. 4, p. 1289-1306, April 2006.

[0102][2] S.G. Mallat, S.G.; Zhang, Z. (1993). "Matching Pursuits with Time-Frequency Dictionaries". IEEE Transactions on Signal Processing. 1993 (12): 3397-3415.

[0103] [3] A. Chambolle, R. A. DeVore, N. Y. Lee, and B. J. Lucier, Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage, IEEE Trans. Image Process., 7 (1998), p. 319-335.

[0104] [4] LASSO (L1 standard): Tibshirani, Robert (1996). "Regression Shrinkage and Selection via the lasso". Journal of the Royal Statistical Society. Series B (methodological). Willy. 58 (1): 267-88.

[0105] [5] M. Aharon, M. Elad and A. Bruckstein, "K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation," in IEEE Transactions on Signal Processing, vol. 54, no. 1 1 , p. 4311-4322, Nov. 2006, doi: 10.1109/TSP.2006.881199.

Claims

CLAIMS Method for reconstructing an image I of an object from a compressed acquisition of measurements, the image I having a number d of pixels, the method comprising the steps of:

- Acquire (201), in a vector m, a number p, strictly less than d, of imaging measurements of the object, by means of an imaging device,

- Determine (202) a first matrix ^AT equal to the product of a projection matrix P defining a set of projection vectors and a dictionary matrix D made up of the components of a decomposition base of the image,

- Find (203) a sparse solution x, of dimension N strictly greater than d, to the underdetermined linear system A ^T x=m in: i. Determining (301 ) a solution x _MC of minimum quadratic norm of said linear system and normalizing the vector x _M c ii. Determining (302) a second matrix whose rows form an orthonormal basis of the complementary subspace orthogonal to the rows of the first matrix ^AT , iii. Searching (303) for a linear combination x of the vector x _M c and of the rows of the second matrix A having a maximum kurtosis.

- Reconstruct (204) the image I by the product of the dictionary matrix D and the vector x, solution of said system. Image reconstruction method according to Claim 1, in which the said linear combination is obtained by searching, for at least one row of the second matrix A , for a linear combination of a pair of components formed by the vector x _M c st of this row , having maximal kurtosis. Image reconstruction method according to Claim 2, in which the said linear combination is obtained by searching for a Givens rotation matrix, defined by a rotation angle θ, for which the product of the said pair of components and of this Givens rotation matrix exhibits a maximum kurtosis for one of its two components or exhibits a maximum difference between the kurtosis of the two components. Image reconstruction method according to Claim 3, in which the search for a linear combination of maximum kurtosis is iterated for several rows of the second matrix A ^T . Image reconstruction method according to any one of the preceding claims, in which the dictionary matrix D consists of a base of elementary signals, for example a base of wavelet functions, or is generated by automatic learning from images d 'learning. Image reconstruction method according to any one of the preceding claims, in which the vector m is acquired by means of an X-ray detector and consisting of the concatenation of several measurement vectors, each measurement vector comprising several measurements carried out by the detector for a different viewing angle of the object. Image reconstruction method according to Claim 6, in which the projection matrix P contains projection vectors respectively defining directions associated with an angle of view. An image reconstruction method according to any preceding claim wherein the I-image is reconstructed in two dimensions or in three dimensions. System for reconstructing an image I of an object from a compressed acquisition of measurements, the image I having a number d of pixels, the system comprising an imaging device for acquiring, in a vector m, a number p, strictly less than d, of object imaging measurements and a processing unit configured for:

Determine a first matrix ^AT equal to the product of a projection matrix P defining a set of projection vectors and a dictionary matrix D consisting of the components of an image decomposition base,

- Find a sparse solution x, of dimension N strictly greater than d, to the underdetermined linear system A ^T x=m in: i. Determining a solution x _M c of minimum quadratic norm of said linear system and normalizing the vector x _M c ii. Determining a second matrix A ^T j_ whose rows form an orthonormal basis of the complementary subspace orthogonal to the rows of the first matrix A ^T , iii. Finding a linear combination x of the vector x _MC and the rows of the second matrix A ^T having maximum kurtosis.

- Reconstruct the image I by the product of the dictionary matrix D and the vector x, solution of the said system. System for reconstructing an image according to claim 9, in which the processing unit is configured to carry out the steps of the method for reconstructing an image according to any one of the preceding claims 1 to 8. System for reconstructing an image An image according to claim 9 wherein the imaging device is an x-ray detector or a magnetic resonance imaging unit. System for reconstructing an image according to one of Claims 9 to 11 comprising a display device for displaying the reconstructed image I.