Abstract:
A method for processing a sequence of a plurality of 2D projection images of an object of interest, acquired with a medical imaging system comprising a source of X-rays adapted to move around the object is provided. The method comprises obtaining 2D projection images of the object of interest according to a plurality of angulations at a first time when the object is without injection of a contrast product; obtaining 2D projection images of the object of interest according to a plurality of angulations at a second time when the object is opacified by injection of a contrast product iteratively treating the projection images by minimizing relative to a sequence of 3D images, wherein the minimizing solution is a set of estimated 3D images.

Description:
BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The present disclosure relates to the field of tomographic reconstruction and more particularly to the reconstruction of a sequence of 3D image(s) describing an object to be imaged by means of a contrast product. 
     2. Description of the Prior Art 
     Tomography allows making images of slices of a region of interest of an object by acquiring projections. 
       FIG. 1  schematically illustrates the acquisition of 2D images of an organ and the tomographic reconstruction of a 3D image of this organ. Reconstruction of a 3D image by tomography consists of emitting X-rays  10  from a source to the organ  12 , the X-rays being emitted according to different angulations lε{1, . . . , L} which define the trajectory T r  of the source (commonly rotation, also known as spin). 
     After having passed through the organ  12 , the X-rays  10  are detected by a detector  13  to form a set of 2D projections. There are as many 2D projections acquired as relevant angulations (or L projections for the trajectory). Acquisition is carried out by a detector (not shown) located opposite the source of X-rays  11 . The detector may be, for example, a digital camera. 
     One application of tomography is the detection and characterisation of a lesion in an organ, for example, stenosis in a vessel of a patient. 
     The acquired 2D projections are used to reconstruct a 3D image of the object. This 3D image is more precisely a 3D map of the coefficients of attenuation to X-rays of the traversed medium. Using this map, the radiological practitioner interprets this image as a function of the differences in contrast observed. 
     An iterative 3D reconstruction process is known. This process is based on a discrete and matricial expression of the problem of tomographic reconstruction. Such a process is carried out 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 set of L projections acquired, R is a projection operator which models the tomographic imaging system and its trajectory of acquisition of L projections p, and f is the 3D image of the object to be reconstructed.
 
     The problem of tomographic reconstruction is determining f by knowing p and R. A known solution to the equation mentioned hereinabove is the resolution of the following criterion: 
                   {             Q   ⁡     (   g   )       =       1   2     ⁢              R   ⁢           ⁢   g     -   p          2   2                     g   *     =         arg   ⁢           ⁢   min     g     ⁢           ⁢     Q   ⁡     (   g   )                         
where ∥ ∥ 2  symbolises the Euclidian standard said L 2 .
 
     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 a few hundreds) and when the set of L angulations covers at least 180° degrees. These conditions are commonly verified when the organ is static during the time necessary to effect a rotation of the imaging system. 
     A problem arises in getting an image of the blood vessels which do not show a marked difference in contrast relative to the surrounding tissue. It is therefore necessary to inject a so-called “contrast” product, for example an iodised product, into the vessels so as to make them more opaque to X-rays and allow them to be displayed as much in 2D projections as in associated 3D reconstructions. The image of the vessels alone is obtained by subtracted angiography where two acquisitions are made: one without contrast product, known as “mask” and noted p M , and a second one, known as “opacified” and noted p O , identical in terms of geometry, but after injection of the contrast product. The order in which these acquisitions are done does not matter. The image of the vessels in the images is obtained by opacified subtraction minus mask. Subtracted tomographic acquisition occurs therefore as the succession of two tomographic acquisitions, masked and opacified, such that an image of the acquisition mask (same geometry) allowing subtraction corresponds to each image of the opacified acquisition. 
     The 3D mask images, noted f M , opacified, noted f O , and subtracted, noted f S , are respectively defined as the solutions to the problems Rf M =p M , Rf O =p O  and Rf S =p O −p M . They are linked by the relationship f S =f O −f M  and by the symbol R common to all three problems, which translates the geometric identity of the spins and consequently the possibility of subtraction of data p O −p M . 
     The contrast product dilutes rapidly in the blood flow. The maximal rotation speed of the tomographic system is generally used to limit the volume of contrast product necessary for opacified acquisition. The angular sampling conditions of the tomography also involve obtaining a few hundreds of images by acquisition, commonly 600. Angiography systems are however limited in speed of rotation, commonly 40°/s, and in image rate, commonly 50 Hz, such that the total number of images acquired, for example 250, is sufficient for the display of strong contrasts such as bones and opacified vessels, but insufficient for weak contrasts, such as soft tissue and perfused tissue. It is said that contrast resolution is limited by sampling. The 3D image subtracted from the contrast product alone (opacified vessels and perfused tissue) suffers under the same limitations of sampling, although it is the result of two acquisitions, simply because these two acquisitions reproduce the same limited sampling enabling subtraction of 2D projections instead of being complementary to enable an increase in contrast resolution. 
     Reference is now made to the general case where the acquisition mask is modelled by the operator R M , and acquisition opacified by the operator R O . The 3D images f M , and f O  are respectively defined by R M f M =p M  and R O f O =p O  and are linked by the relationship f S =f O −f M . As the hypothesis R M =R O , is no longer necessarily done, the subtraction of data can no longer be written; generation of the subtracted 3D image no longer occurs except by subtraction of the opacified 3D image, reconstruction of opacified spin at rapid rotation (defined by R O ), and of the 3D image mask, reconstruction of a spin mask at rapid rotation, or even of any other spin enabling reconstruction of the object of interest without injection of contrast product (defined by R M ). In particular, spin at slower rotation will accumulate the number of projections necessary for detection of weak contrasts of soft tissue in the 3D image mask. This case corresponds however to an increase in the total number of projections acquired and therefore of the dose of X-rays associated with examination, without consequence for contrast resolution of the opacified 3D image which is not improved for perfused tissue. 
     Embodiments of the present invention consist of simultaneously using the two mask and opacified acquisitions to reconstruct the associated opacified and subtracted 3D mask images in order to reduce the number of projections necessary for detection of weak contrasts in the reconstructed 3D images. For this, it is based on temporal interpretation of acquisition with injection of contrast product in considering two temporal points: t O  for the injection and t M  for the mask. It hypothesises that the time changes caused by injection of contrast product are compressible. The time series to be reconstructed is therefore constituted by the vector {right arrow over (f)}={f(t O ),f(t M )}={f O ,f M } for which a time transform {right arrow over (h)}=H t {right arrow over (f)}={h α ,h δ } can be defined where the constituents h α  and h δ  are 3D images and such that one at least of the constituents of H t {right arrow over (f)} is parsimonious. 
     This hypothesis proposes a method of reconstruction for compensating sub-sampling of variations linked to injection, and therefore more particularly improving reconstruction of the subtracted 3D image. A corollary to the hypothesis of compressibility is that the parts not affected by the contrast product are redundant in the acquisitions, and must therefore be determined by the set of available data {right arrow over (p)}={p(t O ),p(t M )}={p O ,p M } according to the relationship R{right arrow over (f)}={right arrow over (p)} with R diagonal operator by block, where each block of the diagonal models spin, or R=diag{R(t M ),R(t O )}=diag{R M ,R O }. So, without making a hypothesis of compressibility on the mask and opacified 3D images themselves, improved reconstruction of their common parts is obtained because of H t . A novel consequence of this analysis is that it is then advantageous to take R(t M )#R(t O ) to make the two acquisitions as complementary as possible. 
     In particular, if R(t M ) samples a circular trajectory with an angular pitch δθ, and if R(t O ) samples a circular trajectory with the same angular pitch δθ, they will be offset by δθ/2 in embodiments of the present invention so that the operator R=diag{R(t O ),R(t M )} corresponds to a sampling of a circular trajectory with an angular pitch of δθ/2. It is evident that angular sampling for the same number of projections is doubled (and therefore the same dose of X-rays) relative to the common case where R(t M )=R(t O ). 
     BRIEF SUMMARY OF THE INVENTION 
     According to an embodiment of the present invention, a method for processing a sequence of a plurality of 2D projection images of an object of interest, acquired with a medical imaging system comprising a source of X-rays adapted to move around the object, is provided. The method comprises obtaining 2D projection images of the object of interest according to a plurality of angulations at a first time when the object is without injection of a contrast product, the 2D projection images being such that R M f M =p M , where p M  are the projection images at the first time, f M  is the 3D image of the object at the first time and R M  is a projection operator which models acquisition by the imaging system at the first time. The method also comprises obtaining 2D projection images of the object of interest according to a plurality of angulations at a second time when the object is opacified by injection of a contrast product, the 2D projection images being such that R O f O =p O , where p O  are the projection images at the second time, f O  is the 3D image of the object at the second time and R O  is a projection operator which models acquisition by the imaging system at the second time. The method further comprises iteratively treating the projection images by minimizing, relative to a sequence of 3D images, the function:
 
 J ({right arrow over ( g )},λ)=λ S ({right arrow over ( g )})+ Q ({right arrow over ( g )}),
 
where
 
     Q({right arrow over (g)}) is a term of fidelity to the measurements acquired; 
     S({right arrow over (g)})=∥Φ xyz H t {right arrow over (g)}∥ 1  is a restriction where H, is a decomposition time such that {right arrow over (h)}=H t {right arrow over (g)} with {right arrow over (h)}={h α ,h δ } combination of the two constituents of {right arrow over (g)} and where Φ xyz {right arrow over (h)}={Φ α h α ,Φ δ h δ } defines two spatial transforms of compression of both constituents of {right arrow over (h)}, restriction, which will improve the determining of variations linked to the injection, supposedly compressible variations, and select parts common to both images g(t O ) and g(t M ) which will be determined at the same time from the projection images at the first time and the projection images at the second time; and 
     λ is a parameter of compressibility. 
     Also, the minimizing solution is a set of 3D images estimated from:
 
{right arrow over (f)}={f O ,f M },
 
where {right arrow over (f)} represents the 3D image of the object at both of the first and second time.
 
     According to another embodiment of the present invention, a medical imaging system is provided. The medical imaging system comprises: an acquisition unit comprising a source of radiation, a sensor for the acquisition of a plurality of 2D projection images of an object; and a processing unit. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Other characteristics and advantages of the invention will emerge from the following description which is purely illustrative and non-limiting and must be considered in conjunction with the attached diagrams, in which: 
         FIG. 1  schematically illustrates the acquisition of 2D images of an organ and the tomographic reconstruction of a 3D image of the organ. 
         FIG. 2  schematically illustrates a medical imaging system according to an embodiment of the present invention; 
         FIG. 3  is a block diagram of a process according to an embodiment of the present invention; and 
         FIG. 4  schematically illustrates steps of a process according to an embodiment of the present invention. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
       FIG. 2  schematically illustrates a medical imaging system  100  for the acquisition of 2D projection images for the reconstruction of a 3D image of an organ. Such a system is used especially for the detection and characterisation of stenoses in vessels. 
     The medical imaging system  100  comprises a support  1  for receiving a patient  12  to be examined, a source  11  designed to emit a beam  10  of X-rays, a detector  13  arranged opposite 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 source  11  of X-rays and the detector  13  are connected for example by means of an arch  15 . 
     The detector  13  can be a semi-conductor image sensor comprising, for example, caesium phosphorous iodide (scintillater) on a matrix of transistor/photodiode of amorphous silicon. Other adequate detectors are: a CCD sensor, a direct digital detector which directly converts X-rays into digital signals. The detector  13  illustrated in  FIG. 2  is planar and defines a flat image surface, while other geometries can naturally be suitable. 
     The control unit  6  controls the acquisition by fixing several parameters such as the dose of radiation to be emitted by the X-ray source and the positioning of the source  11  and of the detector  13 . It is connected to the arch  15  by wire or wireless connection. 
     The control unit  6  can comprise a reading device (not shown), for example a disc reader, a CD-ROM, DVD-ROM reader, or connection ports for reading the instructions of the process for treating an instruction medium (not shown), such as a diskette, a CD-ROM, DVD-ROM, or USB flash drive or more generally by any removable memory medium or even via a network connection. 
     The storage unit  7  is connected to the control unit  6  for recording parameters and acquired images. It is possible to ensure that the storage unit  7  is located inside the control unit  6  or outside the control unit. The storage unit  7  can be formed by a hard drive or SSD, or any other removable and rewritable storage means (USB flash drives, memory cards etc.). The storage unit  7  can be ROM/RAM memory of the control unit  6 , a USB flash drive, a memory card, central server memory. 
     The display unit  8  is connected to the control unit  6  for displaying the acquired images and/or information on the acquisition control parameters. The display unit  8  can be for example a computer screen, a monitor, a flat screen, plasma screen or any type of display device of known type. Such a display unit  8  allows a practitioner to control reconstruction and/or display of the 2D images acquired. 
     The medical imaging system  100  is coupled to a processing system  200 . The processing system  200  comprises a calculation unit  9  and storage unit  10 . 
     The processing system  200  receives the images acquired and stored in the unit memory  4  of the medical imaging system  100  from which it makes a certain number of treatments (see hereinafter), for example reconstruction of a 3D image from 2D images. 
     Transmission of data from the storage unit  4  of the medical imaging system  100  to the calculation unit  9  of the processing system  200  can be done over an internal or external information network or by means of any adequate physical memory medium such as diskettes, CD-ROM, DVD-ROM, external hard drive, USB flash drive, SD card, etc. 
     The calculation unit  9  is for example a computer or computers, a processor or processors, a microcontroller or microcontrollers, a microcomputer or microcomputers, a programmable automaton or automatons, a specific application integrated circuit or circuits, other programmable circuits, or other devices including a computer such as a workstation. 
     By way of variant, the calculator  9  can comprise a reading device (not shown), for example a disc reader, a CD-ROM, DVD-ROM reader, or connection ports for reading the instructions of the process for treating an instruction medium (not shown), such as a diskette, a CD-ROM, DVD-ROM, or USB flash drive or more generally by any removable memory medium or even via a network connection. 
     Also, the processing system comprises a storage unit  14  for storing data generated by the calculation unit  9 . 
     The calculation unit  9  can be connected to the display unit  8  (as in  FIG. 2 ) or else to another display unit (not shown). 
     The image processing method is for example implemented in the treatment unit  200  of the medical imaging system illustrated in  FIG. 2 . The image processing method reconstructs a sequence of two 3D images {right arrow over (f)}={f(t O ),f(t M )} representing an object with two different times: one, noted t O , with injection of contrast, the other, noted t M , without, giving t O &gt;t M  or t O &lt;t M . 
     For each time, there is a set of images of 2D projections p(t n ) with t n  element of {t O , t M }, obtained by means of the medical imaging system moving according to a trajectory about the object in motion (commonly rotation, also known as spin). This trajectory can be different according to the time in question. 
       FIG. 3  shows the acquisition of sets of 2D projection images {right arrow over (p)}={p(t O ),p(t M )} of the object  12  with and without contrast product. In  FIG. 3 , the greyed object corresponds to the object in which a contrast product was injected by means of an injection device  14  of contrast product. 
     The sets of 2D projection images {right arrow over (p)} are, for example, previously acquired and recovered from the storage unit  14  of the processing unit  200  of the medical imaging system  100  and treatment of 2D projection images is carried out in the calculator  9  of the processing unit  200  of the medical imaging system. 
     Each plurality of sequences of 2D projection images p(t n ) is such that R(t n )f(t n )=p(t n ). R=diag{R(t n )} for t n ε{t O ,t M } is the projection operator, diagonal by block, which models sampling conducted by the medical imaging system during acquisition with and without contrast product. It verifies R{right arrow over (f)}={right arrow over (p)}. 
     In current practice, the choice R(t O )=R(t M ) is privileged, for lack of adequate reconstruction. The reconstructed 3D images are affected by sub-sampling streaks. The aim of the invention is to propose a reconstruction method where sub-sampling can be compensated by a priori mathematics which makes it an advantage to have R(t O )≠R(t M ) and thus increases the contrast resolution of all reconstructions. 
     For each time t n ε{t O ,t M } of the object, the functional of the least squares following are defined: 
               Q   ⁡     (     g   ,     t   n       )       =       1   2     ⁢                R   ⁡     (     t   n     )       ⁢   g     -     p   ⁡     (     t   n     )              2   2             
where ∥ ∥ 2  symbolises the Euclidian standard known as L 2  and g a 3D image.
 
     The composite functional associated with {right arrow over (g)}={g(t O ),g(t M )} is defined: 
     
       
         
           
             
               Q 
               ⁡ 
               
                 ( 
                 
                   g 
                   → 
                 
                 ) 
               
             
             = 
             
               
                 
                   Q 
                   ⁡ 
                   
                     ( 
                     
                       
                         q 
                         ⁡ 
                         
                           ( 
                           
                             t 
                             o 
                           
                           ) 
                         
                       
                       , 
                       
                         t 
                         o 
                       
                     
                     ) 
                   
                 
                 + 
                 
                   Q 
                   ⁡ 
                   
                     ( 
                     
                       
                         g 
                         ⁡ 
                         
                           ( 
                           
                             t 
                             M 
                           
                           ) 
                         
                       
                       , 
                       
                         t 
                         M 
                       
                     
                     ) 
                   
                 
               
               = 
               
                 
                   1 
                   2 
                 
                 ⁢ 
                 
                   
                     
                        
                       
                         
                           R 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             g 
                             → 
                           
                         
                         - 
                         
                           p 
                           → 
                         
                       
                        
                     
                     2 
                     2 
                   
                   . 
                 
               
             
           
         
       
     
     The minimizing relative to {right arrow over (g)} of Q({right arrow over (g)}) reconstructs a sequence of 3D images {right arrow over (g)}*={g*(t O ),g*(t M )} of the object with and without contrast product, where the 3D image g*(t O ) is, however, considerably degraded by the reduced number of projection images contained in the sequence of images p(t O ). The image g*(t M ) is also degraded if the number of 2D projections of p(t M ) is the same as for p(t O ), as is the case in common practice which minimizes the total dose of examination. 
     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 multi - slice helical CT” Med Phys.  34 (11):4526-44] contains for example functionals of least weighted squares alternative to the definition of Q(g,t n ) hereinabove. 
     A restriction of parsimony is also defined
 
 S ( {right arrow over (g)} )=∥Φ xyz   H   t   {right arrow over (g)}∥   1  
 
where ∥ ∥ 1  symbolises the standard known as L 1  and H t  is a decomposition time such that {right arrow over (h)}=H t {right arrow over (g)} with {right arrow over (h)}={h α ,h δ } and Φ xyz  is a transform such that Φ xyz {right arrow over (h)}={Φ α h α ,Φ δ h δ } with Φ α  (respectively Φ δ ) transform enabling spatial compression of the 3D image constituent h α  (respectively h δ ) of {right arrow over (h)}. The type of spatial compression transform (transform in wavelets, gradient, identity, . . . ) is defined independently for each constituent according to the hypotheses of compressibility which can be done on each of the constituents. In particular, if one of the constituents is not compressible, the transform which is associated with it is zero (Φ α  or Φ δ  is zero).
 
     The corollary of the compressibility hypothesis linked to injection is a redundancy hypothesis of the information not affected by injection. The aim of decomposition H t  is to combine the mask and opacified 3D images so that the redundant parts, in other words the parts not affected by the contrast product, and therefore common to the two images, are determined from the information available at the same time in p(t M ) and p(t O ). 
     Integration of measurements and of the compressibility hypothesis is done by defining the functional:
 
 J ( {right arrow over (g)},λ )=λ S ( {right arrow over (g)} )+ Q ( {right arrow over (g)} )
 
     A step of the treatment process consists of determining a sequence of images {right arrow over (g)}*(λ)={g*(t O ,λ),g*(t M ,λ)} which minimize the functional J({right arrow over (g)},λ) relative to {right arrow over (g)} for λ fixed and gives an estimation of {right arrow over (f)}={f(t O ), f(t M )}, solution of the problem of perfectly sampled tomographic reconstruction. The estimation {right arrow over (g)}*(λ) is such that the set of restrictions is verified: variations due to injection are determined by the data obtained at the time t O  only, the redundant parts are determined from times t O  and t M , and the missing parts from the compressibility hypothesis. 
     An iterative convex optimization algorithm for minimizing J({right arrow over (g)},λ) is known, for example one of those described or referenced 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 H. “ Fast gradient - based algorithms for constrained total variation image denoising and deblurring problems.” IEEE Trans. Image Process.  18(11):2419-34]. 
     We note A λ  the iteration for minimizing of J({right arrow over (g)},λ) relative to {right arrow over (g)} for λ fixed and A λ   κ [{right arrow over (g)}] the sequence of 3D images results from the application of κ iterations of the algorithm to the sequence of 3D images {right arrow over (g)}. For {right arrow over (g)} 0  fixed, for example a zero sequence, the algorithm is such that 
                 lim     κ   →   ∞       ⁢           ⁢       A   λ   κ     ⁡     (       g   →     0     )         =         arg   ⁢           ⁢   min       g   →       ⁢           ⁢     J   ⁡     (       g   →     ,           ⁢   λ     )               
and in particular such that
 
     
       
         
           
             
               
                 A 
                 λ 
                 κ 
               
               ⁡ 
               
                 ( 
                 
                   
                     g 
                     → 
                   
                   0 
                 
                 ) 
               
             
             ≈ 
             
               
                 
                   arg 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   min 
                 
                 
                   g 
                   → 
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 J 
                 ⁡ 
                 
                   ( 
                   
                     
                       g 
                       → 
                     
                     , 
                     λ 
                   
                   ) 
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               for 
               ⁢ 
               
                   
               
               ⁢ 
               κ 
               ⁢ 
               
                   
               
               ⁢ 
               finished 
               ⁢ 
               
                   
               
               ⁢ 
               and 
               ⁢ 
               
                   
               
               ⁢ 
               sufficiently 
               ⁢ 
               
                   
               
               ⁢ 
               
                 small 
                 . 
               
             
           
         
       
     
     Current practice is fixing a priori an optimal value λ* of λ. This additional knowledge can be circumvented by a process of degressive restriction, that is, by definition of a decreasing sequence of degrees of compressibility Λ={λ 1 , . . . , λ Ξ } such that λ 1 &gt; . . . &gt;λ Ξ ≧0 for which {right arrow over (g)} 0  an estimation {right arrow over (g)}*(Λ, {right arrow over (g)} 0 ) of {right arrow over (f)} solution of the problem of perfectly sampled tomographic reconstruction is determined from the arbitrary sequence, according to: 
     
       
         
           
               
             
               { 
               
                 
                   
                     
                         
                       ⁢ 
                       
                         
                           
                             g 
                             → 
                           
                           0 
                         
                         , 
                         
                           Λ 
                           = 
                           
                             
                               { 
                               
                                 
                                   λ 
                                   1 
                                 
                                 , 
                                 … 
                                 ⁢ 
                                 
                                     
                                 
                                 , 
                                 λΞ 
                               
                               } 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             gives 
                           
                         
                       
                     
                   
                 
                 
                   
                     
                         
                       ⁢ 
                       
                         
                           
                             g 
                             → 
                           
                           ⁡ 
                           
                             ( 
                             
                               λ 
                               1 
                             
                             ) 
                           
                         
                         = 
                         
                           
                             A 
                             
                               λ 
                               1 
                             
                             κ 
                           
                           ⁡ 
                           
                             [ 
                             
                               
                                 g 
                                 → 
                               
                               0 
                             
                             ] 
                           
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           
                               
                             ⁢ 
                             g 
                           
                           → 
                         
                         ⁡ 
                         
                           ( 
                           
                             λ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             ξ 
                           
                           ) 
                         
                       
                       = 
                       
                         
                           
                             A 
                             
                               λ 
                               ξ 
                             
                             κ 
                           
                           ⁡ 
                           
                             [ 
                             
                               
                                 g 
                                 → 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   λ 
                                   
                                     ξ 
                                     - 
                                     1 
                                   
                                 
                                 ) 
                               
                             
                             ] 
                           
                         
                         ⁢ 
                         
                           ∀ 
                           
                             ξ 
                             ∈ 
                             
                               { 
                               
                                 2 
                                 , 
                                 … 
                                 ⁢ 
                                 
                                     
                                 
                                 , 
                                 Ξ 
                               
                               } 
                             
                           
                         
                       
                     
                   
                 
                 
                   
                     
                         
                       ⁢ 
                       
                         
                           
                             
                               g 
                               → 
                             
                             * 
                           
                           ⁡ 
                           
                             ( 
                             
                               Λ 
                               , 
                               
                                 
                                   g 
                                   → 
                                 
                                 0 
                               
                             
                             ) 
                           
                         
                         = 
                         
                           
                             g 
                             → 
                           
                           ⁡ 
                           
                             ( 
                             λΞ 
                             ) 
                           
                         
                       
                     
                   
                 
               
             
           
         
       
     
     To explain the unfolding of the process for obtaining {right arrow over (g)}*(Λ,{right arrow over (g)} 0 ) which is the approximate solution to the problem, E 1 , {right arrow over (g)} 0  and Λ are fixed. E 2 , the sequence {right arrow over (g)}(λ 1 )=A λ     1     κ [{right arrow over (g)} 0 ] is determined which is the sequence minimising the functional J({right arrow over (g)},λ 1 ) relatively to {right arrow over (g)} for λ 1  fixed starting out from {right arrow over (g)} 0 . 
     Next, {right arrow over (g)}(λ ξ )=A λ     ξ     κ [{right arrow over (g)}(λ ξ-1 )] ∀ξξ{2, . . . , Ξ} is determined iteratively E 3 , E 3 ′, E 3 ″ which corresponds to the sequence {right arrow over (g)}(λ ξ ) minimising the functional J({right arrow over (g)},λ ξ ) relative to {right arrow over (g)} for λ ξ  fixed starting out from {right arrow over (g)}(λ ξ-1 ), resulting in E 4 , {right arrow over (g)}*(Λ,{right arrow over (g)} 0 ) the approximate solution to the problem of tomographic reconstruction. The numbers of the steps do not correspond to the diagram which is itself erroneous: the arrow must return to step E3. 
     The number of approximations Ξ and the number of iterations by approximation κ are fixed such that they form the best compromise between the quality of {right arrow over (g)}*(Λ,{right arrow over (g)} 0 ) and the calculation time necessary for generation of {right arrow over (g)}*(Λ,{right arrow over (g)} 0 ) which are both proportional to Ξ and κ. 
     As a variant, the number of iterations κ is advantageously selected depending on the degree of compressibility λ ξ  and inversely proportional to the latter. 
     The sequence of 3D images {right arrow over (g)}*(Λ,{right arrow over (g)} 0 ) therefore represents reconstruction of the object at times t M  and t O , that is, with and without injection of contrast. 
     Transforms and algorithms illustrating the invention as described are explained. The minimizing iteration A λ  of J({right arrow over (g)},λ)=λS({right arrow over (g)})+Q({right arrow over (g)}) can be written as:
 
 A   λ   [{right arrow over (g)}]=proxS   τ=ρλ   [{right arrow over (g)}−ρ∇Q ( {right arrow over (g)} )]
 
where:
 
     ∇Q({right arrow over (g)}) is the gradient of Q({right arrow over (g)}), 
     ρ is a scalar such as ρ∥∇Q({right arrow over (g)})∥&lt;1 ∀{right arrow over (g)}, and 
                 proxS     τ   =     ρ   ⁢           ⁢   λ         ⁡     [     g   →     ]       =         arg   ⁢           ⁢   min       k   →       ⁢     {                k   →     -     g   →            2   2     +     τ   ⁢            S   ⁡     (   k   )            1         }             
is the proximal operator for applying the parsimony restriction.
 
     In addition, if H t  is defined by the Haar transform itself defined by h α =g(t O )+g(t M ) and h δ =g(t O )−g(t M ), it is clear that h, is the subtracted 3D image. If the hypothesis is made that h δ  contains only contrast vessels, h δ  is naturally parsimonious and Φ xyz  can be defined such that Φ xyz H t {right arrow over (g)}={0,h δ } and S({right arrow over (g)})=∥h δ ∥ 1 . This gives proxS τ [{right arrow over (g)}]=H t   −1 {h α ,T τ h δ } where T τ  is a operator of soft thresholding of threshold τ which inserts in between the Haar transform H t  and its inverse H t   −1 . 
     Alternatively, Φ xyz H t {right arrow over (g)}={μ|∇h α |,|∇h δ |} can be defined where ∇ represents the gradient of a 3D image, such that S({right arrow over (g)})=μ|∇h α ∥ 1 +∥∇h δ ∥ 1 . This restriction, known as total variation, is applied to each of the constituents of {right arrow over (h)}. This gives proxS τ [{right arrow over (g)}]=H t   −1 {F τμ h α ,F τ h δ } where F τ and F τμ  are the filters of total parameter variation τ and τμ respectively with μ real positive or zero which modulates the force of total variation following the constituent of {right arrow over (h)}. The restriction therefore returns to filter each constituent of the Haar transform with a filter of total variation of different force. Since the degree of parsimony of h α , an image containing the mask, is supposed to be smaller than the degree of parsimony of h δ , an image of vessels, this results in 0≦μ&lt;1. 
     These two choices are valid irrespectively of the definitions of R(t n ). However, if R(t M ) samples a circular trajectory with an angular pitch δθ, and if R(t O ) samples a circular trajectory with the same angular pitch δθ, it is advantageous with this invention to take R(t M )≠R(t O ) and offsets of δθ/2 so that the operator R=diag{R(t O ),R(t M )} corresponds to a sampling of a circular trajectory with an angular pitch of δθ/2. This choice doubles the angular sampling for the same total number of projections (and therefore the same dose of X-rays) relative to the current common case where R=R(t M )=R(t O ). More generally, if {α M,1 , . . . , α M,N     M   } the set of N M  angulations contained in R(t M ) and {α O,1 , . . . , α O,N     O   } the set of N O  angulations contained in R(t O ) are noted, these sets are selected such that there are two sub-assemblies of N angles {β M,1 , . . . , β M,N }⊂{α M,1 , . . . , α M,N     M   } and {β O,1 , . . . , β O,N }⊂{α O,1 , . . . , α O,N     O   } defining two rotations of angular pitch Δθ offset by an angle 0&lt;φ&lt;Δθ, or: ∀iε{1, . . . , N} β M,i =iΔθ |β M,i −β O,i |=φ. 
     The treatment process of radiological images can be advantageously implemented as a computer program comprising machine instructions for executing the process.