Tomographic reconstruction of a moving object

Abstract

A method is provided for processing a sequence of sets of 2D projection images of a moving object, wherein the sequence of sets of 2D projection images is obtained by a medical imaging system that is in motion along a trajectory. The method comprises determining a sequence of images which minimize a function dependant on a set of 3D images, a term of fidelity of the sequence of images, a function of spatial and temporal compression of the sequence of images, a compressibility parameter, and a sequence of operators for an approximate modelling of motion. The sequence of operators leads to a compression constraint augmented by partial knowledge of the motion and the minimization comprises defining a decreasing sequence of degrees of compressibility for which an estimation is determined from an initial sequence.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

In general, the field of the invention concerns tomographic reconstruction and more particularly the reconstruction of a sequence of 3D image(s) describing a moving object to be imaged.

2. Description of Related Art

Tomography allows images of an object to be obtained in the form of slices of a region of interest, by acquiring projections. FIG. 1 schematically illustrates the acquisition of 2D images of an object such as an organ, and the tomographic reconstruction of a 3D image of this organ. The reconstruction of a 3D image by tomography comprises the emission of X-rays 10 from a source towards the organ 12, the X-rays being emitted at different angles lε{1, . . . , L} which define the trajectory Tr of the source (commonly an angle of rotation). After passing through the organ 12, the X-rays 10 are detected by a detector (not illustrated) to form a set of 2D projections. There are as many 2D projections acquired as angles under consideration (i.e. L acquisitions). Acquisition is initiated by the detector, which is located facing the X-ray source 11, e.g. a digital camera.

One application of tomography is the detection and characterization of a body lesion e.g. stenosis of a patient's vessel. The acquired 2D projections are used to reconstruct a 3D image of the object. This 3D image is more precisely a 3D mapping of the X-ray attenuation coefficients of the exposed medium. It is by means of this mapping that a radiologist practitioner interprets this image in relation to observed differences in contrast.

An iterative 3D reconstruction method is known. This method is based on a discrete matrix expression of the tomographic reconstruction problem. The method is implemented in a processing unit of a medical imaging system. More precisely, the problem is modelled by the following equation:Rf=p where p is the vector of the L acquired projections, R is a projection operator which models the tomographic imaging system and its trajectory for acquisition of the L projections p, and f is the 3D image of the object to be reconstructed. The tomographic reconstruction problem is to determine f having knowledge of p and R.

One known solution to the above-mentioned equation is the solving of the following criterion:

$\hspace{1em}\{\begin{array}{c}Q\left(g\right)=\frac{1}{2}{\mathrm{Rg}-p}_2^2\\ g^{*}=\underset{g}{\arg \min }Q\left(g\right)\end{array}$ where ∥ ∥₂ symbolizes the Euclidian norm called L₂. Minimization of the criterion Q(g) relative to g gives good results (g*≈f) when the set of projections p is such that L is large (typically several hundred) and when the set of L angles covers at least 180°. These conditions are commonly verified when the organ is static during the time needed to perform a rotation of the imaging system.

A problem arises when the organ is in motion and is therefore represented by a series of 3D images {right arrow over (f)}={f(t₁), . . . , f(t_N)} where f(t_n) represents the object at the position referenced t_n. The set of projections p is interpreted as a sequence of sets of projections {right arrow over (p)}={p(t₁), . . . , p(t_N)} where each p(t_n) contains the angles available when the object is in the position referenced t_n. With a single rotation of the system, the reconstruction of f(t_n) is degraded.

if the motion is ignored, and reconstruction is performed by minimizing the criterion Q(g), the solution obtained is:

$g^{*}=\underset{g}{\arg \min }Q\left(g\right)=\sum _{n=1}^{N}f\left(t_n\right)$ It is degraded by object motion which is averaged on the reconstruction. The modelling of the trajectory of the system reduced to the angles of the projections p(t_n) is denoted R(t_n). The 3D image reconstructed from:

$\hspace{1em}\{\begin{array}{c}Q\left(g,t_n\right)=\frac{1}{2}{R\left(t_n\right)g-p\left(t_n\right)}_2^2\\ g^{*}\left(t_n\right)=\underset{g}{\arg \min }Q\left(g,t_n\right)\end{array}$ is degraded since R(t_n) differs from R in that the number of projections of R(t_n), lower than L, is too small and in that the angle coverage is potentially less than 180° C. The sampling condition therefore requires that each p(t_n) contains L projections and covers 180° C. It can then be seen that {right arrow over (p)} must contain N×L projections. Yet, such acquisition cannot be envisaged since it implies multiplying N times the X-ray dose (and optionally the quantity of contrast product) and additionally it lasts N times longer.

One alternative consists of integrating the motion of the object by means of a sequence of operators M={M(t₁), . . . , M(t_N)} which exactly models object motion i.e. such that:∀n,M(t_n)f(t_n)=f(t_ref)where f(t_ref) is the 3D image of the object at a reference position. Reconstruction is then obtained by minimizing the criterion:

$\hspace{1em}\{\begin{array}{c}Q\left(g,M,t_n\right)=\frac{1}{2}{R\left(t_n\right)M^{-1}\left(t_n\right)g-p\left(t_n\right)}_2^2\\ g^{*}\left(t_{\mathrm{ref}}\right)=\underset{g}{\arg \min }\sum _{n=1}^{N}Q\left(g,M,t_n\right)\end{array}$ so that the complete sequence {right arrow over (p)}={p(t₁), . . . , p(t_N)} is used to obtain a single image g*(t_ref) from which each position g*(t_n)=M⁻¹(t_n)g*(t_ref) can be inferred. However, knowledge of M is most often only approximate or only valid in a restricted part of the image. For example, in the document [Blondel C, Malandain G, Vaillant R, Ayache N., “Reconstruction of coronary arteries from a single rotational X-ray projection sequence” IEEE Trans. Med. Imaging 25(5):653-63], the motion of the coronaries is estimated in the projected images and integrated in the reconstruction. However, only the image of the coronaries is obtained, all the other structures being removed from the projective images by pre-processing.

Alternatively, reconstruction is obtained independently for each position t_n by minimizing the criterion:

$\hspace{1em}\{\begin{array}{c}S\left(g\right)={W_{\mathrm{xyz}}g}_1\\ g^{*}\left(t_n,\lambda \right)=\underset{g}{\arg \min }\left\{\lambda S\left(g\right)+Q\left(g,t_n\right)\right\}\end{array}$ where ∥ ∥₁ is the so-called L₁ norm, W_xyz is a spatial transformation enabling the individual compression of the 3D images (e.g. wavelets, gradient, identity (denoted Id . . . ) and λ is a prior fixed scalar. These compressed images have not yet given any convincing results for clinical application. It is on reconstructing solely the coronary vessels, after removing all other information from the projections by pre-processing [Hansis E, Carrol Schäfer D, Dössel O, Grass M., “High Quality 3-D coronary artery imaging on an interventional C-arm X-ray system” Med Phys. 37(4):1601-9] that the criterion W_xyz=Id is known to produce a sparse image of the coronaries i.e. equal to 0 in all the pixels except the vessels whose volume is known to be very restricted (a few percent in the 3D image).

In the document [Chen G H, Tang J, Leng S., “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets” Med. Phys. 35(2): 660-3] it is shown that it is possible to avoid the pre-processing of projection images, provided a strong assumption is made of prior image knowledge denoted f_p, close to the solution. The constraint S(g) is then replaced by αS(g)+(1−α)S(g−f_p) where 0≦α<1. Variants of the cited techniques are compared in the document [Bergner F, Berkus T, Oelhafe M, Kunz P, Pan T, Grimmer R, Ritschl L, Kachelrieβ M, <<An investigation of 4D cone-beam CT algorithms for slowly rotating scanners>> Med. Phys. 37(9): 5044-53].

BRIEF SUMMARY OF THE INVENTION

One objective of the invention is to obtain a 3D reconstruction of a moving object, represented by a series of 3D images {right arrow over (f)}={f(t₁), . . . , f(t_N)} from a sequence of 2D projection images {right arrow over (p)}={p(t₁, . . . , p(t_N)} obtained with a medical imaging system.

The prior art shows that there exist reconstruction algorithms allowing the generation of compression approximations of a 3D image of an object. These reconstructions are obtained by minimizing a criterion of the form

$g^{*}\left(\lambda \right)=\underset{g}{\arg \min }\left\{\lambda S\left(g\right)+Q\left(g\right)\right\}$ where Q(g) models the relationships between the solution and measurement data, S(g) imposes the type of compressibility of the approximate solution, whereas the degree of compressibility is fixed a priori as being proportional to the scalar λ. That in general, the solution g*(λ) is not unique and depends upon the starting point of the minimization procedure. That in general there exists an approximation error g*(λ)−f, itself proportional to λ (it will be said that the solution is biased, however this bias is not uniform on the solution). In addition, that in general it is not possible to have exact prior knowledge of the motion M={M(t₁), . . . , M(t_N)} of the object, but on the other hand it is frequent to have partial prior information on the motion.

A first aspect of the invention is a method for processing a sequence of sets of 2D projection images of a moving object. The motion of which is described by a set of positions referenced by t={t₁, . . . , t_N} wherein the sequence of sets of 2D projection images is obtained by a medical imaging system that is in motion along a trajectory. The method comprises determining a sequence of images which minimize a function dependant on a set of 3D images, a term of fidelity of the sequence of images, a function of spatial and temporal compression of the sequence of images, a compressibility parameter, and a sequence of operators for an approximate modelling of motion. The sequence of operators leads to a compression constraint augmented by partial knowledge of the motion. Minimization comprises defining a decreasing sequence of degrees of compressibility for which, an estimation is determined from an initial sequence.

A second aspect of the invention is a medical imaging system comprising an acquisition unit. The acquisition unit comprises a radiation source and a sensor. The acquisition unit is configured to acquire a sequence of sets of 2D projection images of a moving object while the acquisition unit is in motion along a trajectory, wherein the motion of the moving object is described by a set of positions. The medical imaging system further comprises a computing unit configured to determine a sequence of images which minimizes a function a function dependant on a set of 3D images, a term of fidelity of the sequence of images, a function of spatial and temporal compression of the sequence of images, a compressibility parameter, and a sequence of operators for an approximate modelling of motion.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate one or more embodiments and, together with the description, explain these embodiments. In the drawings:

FIG. 1 schematically illustrates the acquisition of 2D images of an object such as an organ;

FIG. 2 schematically illustrates a medical imaging system according to an embodiment of the invention;

FIG. 3 is a flow diagram of a method according to an embodiment of the invention;

FIG. 4 schematically illustrates steps of a method according to an embodiment of the invention.

DETAILED DESCRIPTION OF THE INVENTION

FIG. 2 schematically illustrates a medical imaging system 100 to acquire 2D projection images for the reconstruction of a 3D image of an object such as an organ. The system may be used for the detection and characterization of stenosis in vessels. The medical imaging system 100 comprises: a support 1 intended to receive a patient (not illustrated) to be examined, in whom it is desired to image an organ 12; a source 11 intended to emit an X-ray beam 10; a detector 13 arranged facing the source 11 and configured to detect the X-rays emitted by the) source 11; a control unit 6; a storage unit 7; and a display unit 8. The X-ray source 11 and the detector 13 are connected by means of a C-arm 15, for example.

The detector 13 may be a semiconductor image sensor, for example comprising caesium iodide phosphor (scintillator) on a transistor/amorphous silicon photodiode array. Other suitable detectors are: a CCD sensor or a direct digital detector which directly converts X-rays into digital signals. The detector 13 illustrated in FIG. 2 is planar and defines a planar image surface, however, other geometries are evidently also suitable.

The control unit 6 is used to control an acquisition by setting several parameters such as the radiation dose to be emitted by the X-ray source and the positioning of the source 11 and detector 13. It is connected to the C-arm 15 via a wire or wireless connection. The control unit 6 may comprise a reader device (not illustrated) e.g. a diskette reader, CD-ROM, DVD-ROM reader or connection ports to read the instructions of the processing method from an instruction medium (not shown) e.g. a diskette, CD-ROM, DVD-ROM, USB key or more generally any removable memory medium or via a network connection.

The storage unit 7 is connected to the control unit 6 to record parameters and acquired images. It is possible to place the storage unit 7 either inside or outside the control unit 6. The storage unit 7 may be formed of a hard disk or SSD, or any other removable, re-write storage means (USB keys, memory cards etc. . . . ). The storage unit 7 may be a ROM/RAM memory of the control unit 6, a USB key, memory card or memory of a central server.

The display unit 8 is connected to the control unit 6 to display acquired images and/or information on acquisition control parameters. For example, the display unit 8 may be a computer screen, a monitor, flat screen, plasma screen or any other type of known display device. The display unit 8 allows the practitioner to control the reconstruction and/or display of acquired 2D images.

The medical imaging system 100 is coupled with a processing system 200. The processing system 200 comprises a computing unit 9 and storage unit 14. The processing system 200 receives acquired images stored in the storage unit 7 of the medical imaging system 100 from which it performs a certain number of processing operations (see below) e.g. a reconstruction of a 3D image from 2D images. The transmission of data from the storage unit 7 of the medical imaging system 100 towards the computing unit 9 of the processing system 200 may take place via an internal or external computer network or using any suitable physical memory medium e.g. diskettes, CD-ROM, DVD-ROM, external hard disk, USB key, SD card etc. . . .

The computing unit 9 is one or more computers for example, or one or more processors, one or more microcontrollers, one or more microcomputers, one or more programmable logic controllers, one or more application-specific integrated circuits, other programmable circuits, or other devices which include a computer such as a workstation. As a variant, the computer 9 may comprise a reader device (not illustrated) for example a diskette reader, CD-ROM or DVD-ROM reader, or connection ports to read the instructions of the processing method from an instruction medium (not illustrated) e.g. a diskette, CD-ROM, DVD-ROM, or USB key or more generally any removable memory medium or via a network connection. In addition, the processing system 200 comprises a storage unit 14 to store data generated by the computing unit 9.

The computing unit 9 may be connected to the display unit 8 (as in FIG. 2) or to another display unit (not illustrated). In addition, the processing unit 200 can be included in the medical imaging system 100, their respective storage units 6,7 in this case being merged.

The method for processing images is implemented, for example, in the processing unit 6 of the medical imaging system 100 illustrated in FIG. 2. The method for processing images allows the reconstruction of a sequence of 3D images {right arrow over (f)}={f(t₁), . . . , f(t_N)} representing a moving object for a set of positions referenced by t={t₁, . . . , t_N}. For each position, there is a set of 2D projection images p(t_n) obtained by means of the medical imaging system moving along a trajectory around the moving object (with the trajectory commonly being an angle of rotation).

FIG. 3 illustrates the acquisition of sets of 2D projection images p(t_n) for different positions t_n of the (moving) object 12. For each instant t_n of the motion, there is a corresponding set of 2D projection images p(t_n). The sets of 2D projection images p(t_n) may be previously acquired, for example, and retrieved from the storage unit 14 of the processing unit 200 or from the storage unit 7 of the medical imaging system 100, and processing of the 2D projection images is performed in the computer 9 of the processing unit 200. Each set of 2D projection images p(t_n) is such that R(t_n)f(t_n)=p(t_n) where R={R(t₁), . . . , R(t_N)} is the sequence of the projection operators which model the sampling made by the medical imaging system for each position t_n of the object while the system is in movement.

For each position t_n of the object, the following least-squares function is defined:

$Q\left(g,t_n\right)=\frac{1}{2}{R\left(t_n\right)g-p\left(t_n\right)}_2^2$ and the associated composite function:

$Q\left(\stackrel{->}{g}\right)=\frac{1}{2}\sum _{n=1}^N{R\left(t_n\right)g\left(t_n\right)-p\left(t_n\right)}_2^2.$ Minimization relative to {right arrow over (g)} of Q({right arrow over (g)}) allows reconstruction of a sequence of 3D images {right arrow over (g)}={g*(t₁), . . . , g*(t_N)} of all the positions of the object, in which each 3D image g*(t_n) is considerably degraded by the reduced number of projection images contained in each sequence of images p(t_n). For example, in the document [Riddell C, Savi A, Gilardi M C, Fazio F, “Frequency weighted least squares reconstruction of truncated transmission SPECT data.” IEEE Trans. Nucl. Sci. 43(4):2292-8] and [Thibault J B, Sauer K D, Bouman C A, Hsieh J., “A three-dimensional statistical approach to improved image quality for multislice helical CT”. Med Phys. 34(11):4526-44] or alternative, weighted least-squares functions can be found to the above equations.

In addition, a function of spatial and temporal compression of the image sequences is defined as follows:S({right arrow over (g)})=∥ℑ_tW_xyz{right arrow over (g)}∥₁ where W_xyz{right arrow over (g)}={W_xyzg(t₁), . . . , W_xyzg(t_N)} and is a spatial transform enabling individual compression of the 3D images (e.g. wavelets, gradient, identity, . . . ) and ℑ_t is a transform (e.g. derived from Fourier analysis) along the axis t followed by the motion.

It is assumed that there is prior knowledge of a sequence of operators M={M(t₁), . . . , M(t_N)} for approximate modelling of motion allowing the construction of an augmented compression constraint such as S(M{right arrow over (f)})≦S({right arrow over (f)}) where M{right arrow over (f)}={M(t₁)f(t₁), . . . , M(t_N)f(t_N)} and {right arrow over (f)} is the solution to the fully sampled problem of tomographic reconstruction.

The integration of measurements and of the compressibility assumption augmented by prior knowledge of the motion, that is approximate, is made by defining the function J({right arrow over (g)},λ)=λS(M{right arrow over (g)})+Q({right arrow over (g)}). One step of the processing method comprises the determination of a sequence of images {right arrow over (g)}*={g*(t₁), . . . , g*(t_N)} which minimizes the function J({right arrow over (g)}, λ) relative to {right arrow over (g)} for fixed λ and gives an approximate solution off {right arrow over (f)}={f(t₁), . . . , f(t_N)}, the solution to the perfectly sampled tomographic reconstruction problem. It is assumed that a convex optimization iterative algorithm is known such as those described for example in [Afonso M V, Bioucas-Dias J M, Figueiredo M A., “Fast image recovery using variable splitting and constrained optimization,” IEEE Trans Image Process. (9):2345-56] and [Beck A, Teboulle M. “Fast gradient-based algorithms far constrained total variation image denoising and deblurring problems.” IEEE Trans Image Process. 18(11):2419-34] for minimization of J({right arrow over (g)}, λ).

The iteration is denoted A_λ allowing the minimization of J({right arrow over (g)}, λ) relative to {right arrow over (g)} for fixed {right arrow over (λ)} and {right arrow over (h)}=A_λ^κ[{right arrow over (g)}] is the sequence of 3D images resulting from the application of κ iteration(s) of the algorithm to the sequence of 3D images {right arrow over (g)}. For fixed {right arrow over (g)}⁰, e.g. a null sequence, the algorithm is such that

$\underset{\kappa ->\infty }{\lim }{A}_l^{\kappa }\left(\stackrel{->0}{g}\right)=\underset{g}{\arg \min }J\left(\stackrel{->}{g},\lambda \right),$ and in particular such that

$A_{\lambda }^{\kappa }\left(\stackrel{->0}{g}\right)\approx \underset{\stackrel{->}{g}}{\arg \min }J\left(\stackrel{->}{g},\lambda \right)$ for κ finite and sufficiently small.

The reconstruction method comprises the definition of a decreasing sequence of degrees of compressibility Λ={λ₁, . . . , λ_E} such as λ₁> . . . >λ_E≧0 for which, from the arbitrary sequence {right arrow over (g)}⁰, an estimation {right arrow over (g)}*(Λ, {right arrow over (g)}⁰) is determined of {right arrow over (f)}, the solution to the perfectly sampled tomographic reconstruction problem, as follows:

$\hspace{1em}\{\begin{array}{c}\stackrel{->0}{g},\Lambda =\left\{{\lambda }_1,\dots ,{\lambda }_{\Xi }\right\}\mathrm{given}\\ \stackrel{->}{g}\left({\lambda }_1\right)=A_{{\lambda }_1}^{\kappa }\left[\stackrel{->0}{g}\right]\\ \stackrel{->}{g}\left({\lambda }_{\xi }\right)=A_{{\lambda }_{\xi }}^{\kappa }\left[\stackrel{->}{g}\left({\lambda }_{\xi -1}\right)\right]\forall \xi \in \left\{2,\dots ,\Xi \right\}\\ \stackrel{->*}{g}\left(\Lambda ,\stackrel{->0}{g}\right)=\stackrel{->}{g}\left({\lambda }_{\Xi }\right)\end{array}$ With each iteration, the algorithm allows determination of the sequence of 3D images {right arrow over (g)} which, for a given value of λ_ξ, minimizes the function or J({right arrow over (g)},λ_ξ).

As a variant, the number of iterations κ is preferably chosen to depend on the degree of compressibility λ_ξ and is inversely proportional thereto.

FIG. 4, diagrams the method to obtain {right arrow over (g)}*(Λ, {right arrow over (g)}⁰), which is the approximate solution to the problem. In a first step E₁, {right arrow over (g)}⁰ and Λ are fixed. In a second step E₂, a sequence {right arrow over (g)}(λ₁)=A_λ1^κ[{right arrow over (g)}⁰] is determined, which is the sequence that minimizes the function J({right arrow over (g)}, λ₁) relative to {right arrow over (g)} for a fixed λ₁ starting from{right arrow over (g)}⁰. Next, iteratively, in another steps E₃, E₃′, E₃″ we determine {right arrow over (g)}(λ_ξ)=A_λξ^κ[{right arrow over (g)}(λ_ξ−1)]∀ξε{2, . . . , Ξ} which corresponds to the sequence {right arrow over (g)}(λ_ξ) that minimizes the function J({right arrow over (g)}, λ_ξ) relative to {right arrow over (g)} for a fixed λ_ξ starting from {right arrow over (g)}(λ_ξ−1). In step E₄, {right arrow over (g)}*(Λ, {right arrow over (g)}⁰) is obtained, which is the approximate solution to the problem of tomographic reconstruction. The number of approximations Ξ and the number of iterations per approximation κ are fixed such that they form the best compromise between the quality of {right arrow over (g)}*(Λ,{right arrow over (g)}⁰) and the computing time needed to generate {right arrow over (g)}*(Λ, {right arrow over (g)}⁰), which is proportional to both Ξ and κ. The sequence of 3D images {right arrow over (g)}*(Λ,{right arrow over (g)}⁰) therefore represents the reconstruction of the moving object.

The method for processing radiological images can advantageously be implemented in the form of a computer programme comprising machine instructions for applying the method.

Claims

1. A method for processing a sequence of sets of 2D projection images of a moving object, the motion of which is described by a set of positions referenced by t={t1, . . . , tN}, wherein the sequence of sets of 2D projection images is obtained by a medical imaging system that is in motion along a trajectory, wherein each set of 2D projection images is defined by the following equation: R(tn)f(tn)=p(tn),

2. The method according to claim 1 wherein prior knowledge of motion is such that M(tn) is the identity.

3. The method according to claim 1, wherein Q({right arrow over (g)}) is a criterion of ordinary, of weighted or of filtered least-squares.

4. The method according to claim 1, wherein S({right arrow over (h)})=∥ℑtWxyz{right arrow over (h)}∥1 where Wxyz{right arrow over (h)} a spatial transform allowing the individual compression of each 3D image of the sequence of 3D images and where Wxyz{right arrow over (h)}={Wxyzh(t1), . . . , Wxyzh(tN)}, t is an axis, ℑ, is a transform derived from a, Fourier analysis along the axis following the motion of the moving object and ∥ ∥1 symbolizes the L1 norm.

5. The method according to claim, 1 wherein the initial sequence is a null sequence.

6. The method according to claim 1, wherein the sequence of sets of 2D projection images are previously acquired images stored in a storage unit of the medical imaging system.

7. A medical imaging system comprising: an acquisition unit comprising a radiation source and a sensor, the acquisition unit configured to acquire a sequence of sets of 2D projection images {right arrow over (p)}={p(t1), . . . , p(tN)} of a moving object while the acquisition unit is in motion along a trajectory, wherein the motion of the moving object is described by a set of positions referenced by t={t1, . . . , tN }, wherein each set of 2D projection images is defined by the following equation: R(tn)f(tn)=p(tn),