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.
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:
where ∥ ∥2 symbolises the Euclidian standard said L2.
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 pM, and a second one, known as “opacified” and noted pO, 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 fM, opacified, noted fO, and subtracted, noted fS, are respectively defined as the solutions to the problems RfM=pM, RfO=pO and RfS=pO−pM. They are linked by the relationship fS=fO−fM 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 pO−pM.
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 RM, and acquisition opacified by the operator RO. The 3D images fM, and fO are respectively defined by RMfM=pM and ROfO=pO and are linked by the relationship fS=fO−fM. As the hypothesis RM=RO, 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 RO), 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 RM). 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: tO for the injection and tM 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(tO),f(tM)}={fO,fM} for which a time transform {right arrow over (h)}=Ht{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 Ht{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(tO),p(tM)}={pO,pM} 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(tM),R(tO)}=diag{RM,RO}. 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 Ht. A novel consequence of this analysis is that it is then advantageous to take R(tM)#R(tO) to make the two acquisitions as complementary as possible.
In particular, if R(tM) samples a circular trajectory with an angular pitch δθ, and if R(tO) 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(tO),R(tM)} 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(tM)=R(tO).
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 RMfM=pM, where pM are the projection images at the first time, fM is the 3D image of the object at the first time and RM 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 ROfO=pO, where pO are the projection images at the second time, fO is the 3D image of the object at the second time and RO 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)})=∥ΦxyzHt{right arrow over (g)}∥1 is a restriction where H, is a decomposition time such that {right arrow over (h)}=Ht{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(tO) and g(tM) 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)}={fO,fM},
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.
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:
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
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
The image processing method is for example implemented in the treatment unit 200 of the medical imaging system illustrated in
For each time, there is a set of images of 2D projections p(tn) with tn element of {tO, tM}, 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.
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(tn) is such that R(tn)f(tn)=p(tn). R=diag{R(tn)} for tnε{tO,tM} 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(tO)=R(tM) 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(tO)≠R(tM) and thus increases the contrast resolution of all reconstructions.
For each time tnε{tO,tM} of the object, the functional of the least squares following are defined:
where ∥ ∥2 symbolises the Euclidian standard known as L2 and g a 3D image.
The composite functional associated with {right arrow over (g)}={g(tO),g(tM)} is defined:
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*(tO),g*(tM)} of the object with and without contrast product, where the 3D image g*(tO) is, however, considerably degraded by the reduced number of projection images contained in the sequence of images p(tO). The image g*(tM) is also degraded if the number of 2D projections of p(tM) is the same as for p(tO), 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,tn) hereinabove.
A restriction of parsimony is also defined
S({right arrow over (g)})=∥ΦxyzHt{right arrow over (g)}∥1
where ∥ ∥1 symbolises the standard known as L1 and Ht is a decomposition time such that {right arrow over (h)}=Ht{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 Ht 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(tM) and p(tO).
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*(tO,λ),g*(tM,λ)} 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(tO), f(tM)}, 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 tO only, the redundant parts are determined from times tO and tM, 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
and in particular such that
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> . . . >λΞ≧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:
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, E1, {right arrow over (g)}0 and Λ are fixed. E2, the sequence {right arrow over (g)}(λ1)=Aλ
Next, {right arrow over (g)}(λξ)=Aλ
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 tM and tO, 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)})∥<1 ∀{right arrow over (g)}, and
is the proximal operator for applying the parsimony restriction.
In addition, if Ht is defined by the Haar transform itself defined by hα=g(tO)+g(tM) and hδ=g(tO)−g(tM), 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 ΦxyzHt{right arrow over (g)}={0,hδ} and S({right arrow over (g)})=∥hδ∥1. This gives proxSτ[{right arrow over (g)}]=Ht−1{hα,Tτhδ} where Tτ is a operator of soft thresholding of threshold τ which inserts in between the Haar transform Ht and its inverse Ht−1.
Alternatively, ΦxyzHt{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)}]=Ht−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≦μ<1.
These two choices are valid irrespectively of the definitions of R(tn). However, if R(tM) samples a circular trajectory with an angular pitch δθ, and if R(tO) samples a circular trajectory with the same angular pitch δθ, it is advantageous with this invention to take R(tM)≠R(tO) and offsets of δθ/2 so that the operator R=diag{R(tO),R(tM)} 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(tM)=R(tO). More generally, if {αM,1, . . . , αM,N
The treatment process of radiological images can be advantageously implemented as a computer program comprising machine instructions for executing the process.
Number | Date | Country | Kind |
---|---|---|---|
10 61185 | Dec 2010 | FR | national |
Number | Name | Date | Kind |
---|---|---|---|
20090074277 | Deinzer et al. | Mar 2009 | A1 |
20090161932 | Chen | Jun 2009 | A1 |
20100166267 | Zhang et al. | Jul 2010 | A1 |
Number | Date | Country |
---|---|---|
101553171 | Oct 2009 | CN |
101558293 | Oct 2009 | CN |
102008007830 | Aug 2009 | DE |
2010233762 | Oct 2010 | JP |
Entry |
---|
French Search Report dated Apr. 12, 2011 which was issued in connection with the French Application No. 1061185 which was filed on Dec. 23, 2010. |
Chen Guang-Hong et. al.: “Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets”, Medical Physics, AIP, Melville, NY, US, vol. 35, No. 2, Jan. 28, 2008, pp. 660-663, XP012115917, ISSN: 0094-2450, DOI: DOI: 10.1118/1.2836423. |
Mark Supanich et. al.: “Radiation dose reduction in time-resolved CT angiography using highly constrained back projection reconstruction; Radiation dose reduction in time-resolved CTA using HYPR reconstruction” , Physics in Medicine and Biology, Taylor and Francis Ltd. London, GB, vol. 54, No. 14, Jul. 21, 2009, pp. 4575-4593, XP020158912, ISSN: 0031-9155. |
Search Report and Written Opinion from corresponding EP Application No. 11193790 dated Apr. 19, 2012. |
Chen Guang-Hong et al., “Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets”, Medical Physics, vol. 35, No. 2, pp. 660-663, Jan. 28, 2008. |
Mark Supanich et al., “Radiation dose reduction in time-resolved CT angiography using highly constrained back projection reconstruction; Radiation dose reduction in time-resolved CTA using HYPR reconstruction”, Physics in Medicine and Biology, vol. 54, No. 14, pp. 4575-4593, Jul. 21, 2009. |
Riddell et al., “Frequency Weighted Least Squares Reconstruction of Truncated Transmission SPECT Data”, IEEE Transactions on Nuclear Science, vol. 43, Issue No. 4, pp. 2292-2298, Aug. 1996. |
Afonso et al., “Fast Image Recovery Using Variable Splitting and Constrained Optimization”, IEEE Transactions on Image Processing, vol. No. 19, Issue No. 9, pp. 2345-2356, Sep. 2010. |
Thibault et al., “A Three-Dimensional Statistical Approach to Improved Image Quality for Multislice Helical CT”, Medical Physics, vol. No. 34, Issue No. 11, pp. 4526-4544, Nov. 2007. |
Beck et al., “Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems”, IEEE Transactions on Image Processing, vol. No. 18, Issue No. 11, pp. 2419-2434, Nov. 2009. |
Unofficial translation of Chinese Search Report from corresponding CN Patent Application No. 2011104632488 dated Sep. 24, 2014. |
Number | Date | Country | |
---|---|---|---|
20120163532 A1 | Jun 2012 | US |