Patent Application
20070127801
The invention relates to computer limited angle computer tomography. Especially advantageously the method according invention is used for X-ray or gamma ray tomography.
Unlike conventional X-ray imaging, 3D (three-dimensional) imaging enables measuring distances and defining exact structures of the objects. For example in dental implantology, it is critical for the firm attachment of the implant to use optimal screw size and therefore measure optimal depth and angle for the screw hole If the bone implant screw is too short, the implant will be loose. On another hand, if the hole is drilled too deep, the mandibular nerve could be damaged. Other example of the benefit in 3D imaging is avoiding the superposition problem (i.e. problem caused by overlapping issues) in diagnostic imaging.
So far the medical 3D imaging has been associated to CT (computed tomography) algorithms and systems. The basic philosophy of CT imaging is to measure the attenuation in each volume element (voxel) cross the object of interest. This approach requires specially designed system with exact mechanical accuracy and over-sampling of the issue of interest, which means in practice high dose and high device costs comparing to 2D (two-dimensional) imaging. Therefore, the use of CT imaging is now-a-days limited to serious diseases and it is typically utilized only in big hospital districts which eventually limits the smooth patient workflow.
Lately there has been increasing interest in dedicated dental cone beam CT-based systems (CBCT). Since the dental CBCT systems are based on same algorithms and philosophy as non-dedicated CTs, these CBCT systems has also the same drawbacks than non-dedicated CTs. Specially, as in dental CBCT system the patient is sitting (instead of laying) during the exposure, there is even more patient motion artifacts in the images than in non-dedicated CT. It is well known that these kinds of artifacts cannot be handled in current CT algorithms. Moreover, despite the fact that the dedicated dental CBCTs are less expensive than conventional CTs, the total cost of the dental CBCT device is typically more than ten times the cost of standard digital panoramic device.
The object of the invention is to improve an algorithm, which could calculate a 3D model from small number of images and could handle incomplete imaging geometry and intensity information in reasonable calculating time. With that kind of construction, the 3D images could be taken by using any existing digital X-ray system, including dental panoramic and intra oral systems or mammography system. Furthermore, since there would be no need for high dose or mechanical accuracy, 3D imaging could be expanded to any medical study and it would fit smoothly into clinical workflow.
Also the invention aims to enable lower cost of the hardware for calculating the tomography. Present limited angle tomography is either very slow or it need expensive dedicated hardware, and still with for example large memory display controller hardware for calculations the calculating time is too slow for smooth clinical workflow. The aim is to use standard affordable computer systems and still get the results much faster than before.
Biopsy in mammography is another especially advantageous use in addition to all the present uses of computer tomography, as the method enables smaller doses with acceptable image quality and much faster calculation than in present methods.
State of Art, Arithmetic Reconstruction Technique (ART)
One of the most used back projection algorithm in limited angle cases is called Arithmetic Reconstruction Technique (ART). ART is a common name for set of iterative reconstructions where system is modeled with number of linear equations. The ultimate purpose of the ART calculation is to solve these equations.
In limited angle problem a relatively small number (typically from 7 to 11) X-ray images are taken from different angles from the same object. Since there are fewer measurements than variables and contrast and spatial limited accuracy, it is impossible to define the object non-ambiguously. Therefore, instead of giving exact solution, the function of the ART is to minimize the error between the measurement (projection images) and projected volume. This kind of under-determined system is commonly known as an ill-posed system.
The ART algorithm is always a heavy iterative algorithm and the known ART methods are too slow to calculate with present computers. Even the method itself would give satisfactory result; the calculating time is limiting clinical use and there have been no commercial devices available yet.
In known ART algorithm we combine the volume xεR3 and the measurement mεR2 by using geometry matrix A:R3→R2 and volume. If we ignore the noise, we can say
m==Ax →(2.1)
where m is the measured values in matrix, A presents the imaging geometry matrix and x is the 3D presentation of the object. Typically the matrix A is not formed as a matrix, because it is enormously large and sparse. The A is usually presented as a set of rules that are applied for each element when needed. The purpose of reconstruction is to define the x when the m is known. Hence the A−1 does not exist; we need a process to minimize the error between the measurement m and projected volume Ax. However, in X-ray imaging we only know the sum cross the attenuation and therefore, it is not appropriate to use directly the equation ε=Ax−m for likelihood error. Instead of that, we multiply the equation 2.1 by AT:R2→R3 from left-hand side.
ATm=ATAx (2.2)
The operator AT back-projects the measurements into the volume and x is multiplied by operator AT A: R3→R3. Here AT A defines the freedom for x cross the beam (i.e. only the sum cross the beam has to match the measurement not single values). Now we can refine the likelihood error more sensible way
ε=ATAx−ATm (2.3)
If we update the guess x by subtracting the attenuated error from the guess on each iteration round, we get an iterative process which minimizes the error ε, i.e.
xi+1=xi−λεi (2.4)
where i is the iteration round and λ is the relaxation parameter (0<λ<<1). Since it is more effective to do the iteration by using one projection image per iteration round, we rather define
xi+1=xi−λAkT(Akxi−mk) (2.5)
This is one of the most known form of ART. Here k is the projection number (i=Nk, NεN+). The number N is the number of iterations.
Problem of ART Algorithm and its Efficiency in Computer Iteration
The problem of ART algorithms is long computing time. The main reason for that is the non-optimized PC architecture for the sparse operations like Akxi. Technically speaking, the bottleneck in PC architecture is the bus between the main memory (RAM) and the microprocessor. Several architectural fixes has been done to solve that problem. One standard solution for this problem is method called cache, which is a buffer between the memory and the microprocessor. The purpose of the cache is to make sophisticated guess for the next access location on the memory. This is extremely powerful when same or next memory location is accessed several times. However, since matrix A is a sparse matrix and not affine, the cache fails in most of the cases. The effect of the cache is defined by hit rate factor. In typical ART calculation, the cache hit rate is low, which leads to longer memory access times and finally, to time consuming iteration rounds. This phenomenon does not depend whether the volume is updated voxel-by-voxel or using back projection cross the beam. Also there is no cheap computer architecture for fast and random memory access to large memory, the access time to very large memory is always considerably longer than the processors handling time for the data. So there is not expectable in near future hardware solution for the problem.
Method According to the Invention in Affine Imaging Geometry
The invention solves the above mentioned problem with low computer cache hit rate and it even enables distributed calculation. The invention uses an iterative algorithm, which minimizes the likelihood in the frequency domain instead of spatial domain as in prior art. Typically the frequency behavior of the system is defined by calculating Fourier transform of systems PSF (Point Spread Function). In affine case (i.e. all x-ray beams are parallel) AT A can be considered as a convolution and consequently the PSF of the AT A can be defined. Now we can transform equation 2.5 to frequency domain
Xi+1=Xi−λ(Fkxi−Mk) (4.1)
where Fk=FFT(AkTAk), Mk=FFT(AkT mk) and Xi=FFT(xi). Despite the fact that we need to calculate the FFT, as well as AT m, before the iteration is this method significantly faster than using equation 2.5, because the voxel-wise multiplying is applied (marked with {circle around (×)}) and therefore corresponding voxels are on same location and the cache would speed-up the process dramatically. However, in non-affine case the PSF varies in function of space and therefore it is impossible to define the exact convolution kernel for AT A. Therefore, this type of approach is limited only to affine cases.
Method According to the Invention for Arbitrary Imaging Geometry
As we show earlier the straight forward Fourier conversion is limited only in affine cases. Therefore we have to develop another approach for arbitrary geometry. Projection mεR2 from the object xεR3 can be defined as follows:
m(v′)=∫Lx(v) (5.1)
where x(v) is (logarithm of) attenuation in point vεR3 and v′εR2 is the corresponding projection point. Fourier transform of the measurement m is
M(ω1)=∫A∫Lx(v)exp(−i2πωTv′) (5.2)
where A is the area in detector m. On the another hand, the Fourier transform of the object is
X(ω2)=∫Vx(v)exp(−i2π(ω2Tv)) (5.3)
where V is the volume. Notice that ω1εR2 and ω2 εR3.
To compare M and X in equations 5.2 and 5.3 we define two new operators φ:R2→R3 and Λ: R3→R3 so that equivalency φ(M(ω1))=Λ(X(ω2)) exists like shown in
φ is an operator which re-maps the projection image to volume space in frequency domain. Therefore, by using equation 2.2, we elaborate the operator φ
Φ(mk)=FFT(ATmk) (5.4)
This means that the operator φ defines the correlation between the volume and the measurement in frequency domain based on simple back-projection. Despite AT m fulfills the mathematical equations, it typically (almost always) gives the wrong solution, which can be seen as a blurring effect in reconstruction. However, this blurring can be minimized (or at least, controlled) by the operator Λ.
The purpose for the operator Λ is to cut-off all the frequencies outside region defined by φ. In addition, Λ can be considered as a probability matrix; when the corresponding frequency in φ(mk) is known, the value is set to 1, otherwise it is 0. Furthermore, if the geometry is not well-known the values can be set also between 1 and zero. Notice that Λi has the same dimension than Xi and, unlike φ, Λ itself is an irreversible operator.
Now we can define the ART equation respective to equation 2.5 in Fourier domain.
Xi+1=Xi−λΛk(Xi−Φ(mk)) (5.5)
Benefits against the spatial ART (ref equation 2.5) are the voxel-wise multiplying and the information accumulation. For example, in cone beam geometry the well-known frequencies are located on a surface of a ball. The radius of the ball corresponds the SID (source to image distance) and the center point corresponds with the X-ray source point respectively. In fact, this rule is valid in any one-shot imaging geometry.
Test Results
The
The volumes were reconstructed by using different Λ values while other parameters were fixed (including initial guess). Results can be seen on
Therefore the frequency based ART according the invention works well and is remarkably faster than known ART based on spatial domain iterations and the quality of tomography image is not deteriorated compared to the known ART methods.
The present invention has been described in terms of the preferred embodiment, and it is recognized that equivalents, alternative, and modifications, aside from those expressly stated, are possible and within the scope of the appending claims. For example the instead of the Fast Fourier Transform may be used other transforms, like z-transform, wavelet-transform, or any imaginary or real transform, that enables the iteration to be calculated using inner product or sequential computation.
