Method for image reconstruction in a computed tomography apparatus

Abstract

An image reconstruction method for a computer tomograph produces arbitrary eccentric segments of images. The method is a modified linogram reconstruction, more precisely a combination of linogram and gridding methods, by means of which it is possible to reconstruct arbitrary eccentric segments directly, by producing points directly on a Cartesian grip, without the necessity of transformations or interpolations.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention is directed to a method for generating a displayable image from data obtained by irradiation of an examination subject in a computed tomography apparatus, and in particular to a method for image reconstruction suitable for use in such a computed tomography apparatus.

2. Description of the Prior Art

In computed tomography, CT images can be reconstructed by the linogram method. The linogram method was introduced by P. Edholm in 1987, see P. R. Edholm, G. T. Herman, "Linograms in Image Reconstruction from Projections,"IEEE Trans. on Med. Imag., vol. MI-6, No. 4, 1987, pp. 301-307. It is a special Fourier reconstruction method that functions entirely without interpolation in the frequency space. A disadvantage of the usual Fourier reconstruction method is that after the one-dimensional fast Fourier transformation of the measured projections f.sub.lk, discrete spectrum values are defined on a polar grid. The transition to the Cartesian grid in the frequency space takes place by means of a two-dimensional interpolation, which leads to flawed images if insufficient care is taken. The essential idea of the linogram method is to organize discrete spectrum values on a "quasi-Cartesian" grid from the outset through suitable measures, so that interpolation in the frequency space is no longer necessary.

In the known linogram method, a spectrum is always generated that corresponds to the image in the overall measurement field with the diameter D.sub.M. As in the direct Fourier transformation methods, this results in that a very large number of points must be reckoned with in the frequency space with matrices. In the original linogram method, spectral values with a fixed grid ##EQU1## result in the .rho..sub.x direction in the angle region [-45.degree., 45.degree.] in the frequency space, while spectral values with a fixed grid ##EQU2## result in the .rho..sub.y direction in the angle region [45.degree., 135.degree.] (.rho..sub.x and .rho..sub.y are the Cartesian frequency variables). If there are N.sub.p =2048 projections in the angle region [-45.degree., 135.degree.] with N=1024 measurement values per projection in the sampling grid ##EQU3## this means that the two sub-matrices for [-45.degree., 45.degree.] or [45.degree., 135.degree.] in the frequency space each have ##EQU4##

If the method is to be of practical importance, it must be possible to reconstruct arbitrary eccentric image segments. For this purpose, in the known linogram method chirp-z transformations are undertaken on the desired image segment D.sub.B .multidot.D.sub.B with N.sub.Pix .multidot.N.sub.pix image elements in the grid ##EQU5## on the basis of the two sub-matrices, in both the .rho..sub.x and the .rho..sub.y directions. This step is so expensive in terms of processing outlay that it far outweighs the advantage of the lack of interpolation in the frequency space.

SUMMARY OF THE INVENTION

An object of the present invention is to provide a Fourier reconstruction method for image reconstruction for computed tomograph with which arbitrary eccentric image segments can be directly reconstructed.

This object is inventively solved in accordance with the principles of the present invention in a method which is a modified linogram reconstruction method, with which the problem of the representation of an eccentric segment of a CT image is solved with a greatly reduced expense. It is a combination of linogram and gridding methods (on gridding methods, see J. D. O'Sullivan, "A Fast Sinc Function Gridding Algorithm for Fourier Inversion in Computer Tomography,"IEEE Trans. on Med. Image., vol. MI-4, no. 4, 1985, pp. 200-207), with which it is possible to reconstruct arbitrary eccentric image segments directly. The matrix that thereby results in the frequency space has, instead of ##EQU6## or ##EQU7## as in the original linogram method, the frequency grid ##EQU8## (D.sub.B is the side length of the desired image segment). In addition, the points lie on a Cartesian grid, not on the "quasi-Cartesian" grid of the linogram method. If the image consists of N.sub.Pix .multidot.N.sub.Pix =512.multidot.512 image elements, then the Cartesian frequency matrix has 2N.sub.Pix .multidot.2N.sub.Pix =1024.multidot.1024 points, independent of the choice of D.sub.B. In contrast to the original, known gridding method, in the inventive method the gridding itself is only one-dimensional. The back-transformation of the spectrum into the spatial domain can be carried out with fast Fourier transformations; chirp-z transformations, as in the original linogram method, are not necessary.

DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic illustration of a computed tomography apparatus wherein data are obtained, and an image reconstruction is undertaken, in accordance with the principles of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

The computed tomography apparatus shown in the drawing has a measurement unit consisting of an X-ray source 1, which emits a fan-shaped X-ray beam 2, and a radiation receiver 3, consisting of a row of individual detectors, e.g. 1024 individual detectors. The focus is designated 11. The patient 4 to be examined lies on a patient bed 5. For the scanning the patient 4, the measurement unit formed by the source 1 and the receiver mounted on a ring 14 3 is rotated through 360.degree. around a measurement field 9 in which the patient 4 lies by a rotation drive 15. The axis of rotation is designated 10. The X-ray source 1, fed by a high voltage 6, is operated in a pulsed manner or with continuous radiation. At predetermined rotational angle positions of the measurement unit, sets of data are produced that are supplied from the receiver 3 to a computer 7, which calculates the attenuation coefficients of predetermined image points from the registered data sets and graphically reproduces the coefficients on a viewing apparatus 8. An image of the transilluminated layer of the patient accordingly appears on the viewing apparatus 8.

The starting point for the reconstruction is a parallel data set with N.sub.P projections in the angle area [-45.degree., 135.degree.] and N channels per projection in the scanning grid a. If a CT scanner measures data in fan geometry, the parallel data are obtained by rebinning. If f.sub.LK is the measurement value at the projection l in the channel k, then:

.theta..sub.l =-45.degree.+l.DELTA..theta..sub.l l=0, 1, 2, . . . , N.sub.P -1 (1) p.sub.k =(k-0.5+a.sub.m)a k=. . . , -1, 0, 1, . . . (2) wherein a.sub.m is the alignment factor.

Furthermore, .DELTA. tan .theta. or .DELTA. cot .theta. is constant: ##EQU9## The only image that can be calculated from these measurement values without additional assumptions is defined as follows: ##EQU10## wherein M.sub.A (.rho.) is the modulation transfer function of the reconstruction.

Furthermore, one can set ##EQU11## L.sub.C and L.sub.S do not form the basis for restrictions as in a convolution back-projection reconstruction, where L.sub.C is the CT convolution kernel and L.sub.S is the interpolation function necessary in the back projection; however, L.sub.S is a smoothing function with little extension in the spatial domain.

The basic idea of the inventive method is to carry out the necessary convolution of the measurement values f.sub.lk with the CT convolution Kernel in the frequency space as a multiplication, and then to center the image in the desired reconstruction center -r.sub.z as the new origin by means of multiplication by a phase factor. After the one-dimensional back-transformation of the measurement values into the spatial domain, the subsequently convoluted, centered projections can be limited to the desired image field D.sub.B. This is followed by one-dimensional chirp-z transformations to the desired frequency grid ##EQU12## for the angle area [-45.degree., 45.degree.], or ##EQU13## for the angle area [45.degree., 135.degree.]. This is theoretically grounded below.

If the reconstruction center is chosen at -r.sub.z, one first forms B.sub.0 (r-r.sub.z) as: ##EQU14## Consider the convoluted, "centered" projections ##EQU15## The function f.sub.c (p.sub.q, .theta..sub.l) can be periodically repeated without degrading the image in an area with diameter D.sub.M (now however centered around -r.sub.z), if for the period length w, w.gtoreq.2D.sub.M =2Na holds, and if L.sub.c (p) for .vertline.p.vertline. >D.sub.M is set to 0, i.e. ##EQU16## If the image is to be reconstructed in a square centered around -r.sub.z with side length D.sub.B .multidot.D.sub.B, with D.sub.B .ltoreq.D.sub.M, then the function f'.sub.c (p.sub.q, .theta..sub.l) can be limited to this square through multiplication by a rectangular window R.sub.l (p.sub.q), since values outside it do not contribute to the image. A rectangular window R.sub.l (p.sub.q), with extension v.sub.l dependent on the angle of projection .theta..sub.l, is chosen. This is possible as long as the following holds: ##EQU17## In particular, one can set

v.sub.l =2D.sub.B cos .theta..sub.l -45.degree..ltoreq..theta..sub.l .ltoreq.45.degree. v.sub.l =2D.sub.B sin .theta..sub.l -45.degree..ltoreq..theta..sub.l .ltoreq.135.degree. (13) The function f'.sub.c (p.sub.q, .theta..sub.l), limited to v.sub.l with the rectangular window R.sub.l (p.sub.q) and convoluted with L.sub.S, can be periodically repeated with period length v.sub.l without influencing the image in the square D.sub.B .multidot.D.sub.B. This results in ##EQU18## The Fourier-transform B.sub.LG (.rho..sub.x, .rho..sub.y) of B.sub.LG (r-r.sub.z) is as follows: ##EQU19##

As in the linogram method, B.sub.LG (r-r.sub.z) is decomposed into a partial image B.sub.LG,1 (r-r.sub.z) with -45.degree..ltoreq..theta..sub.l .ltoreq.45.degree. and into a partial image B.sub.LG,w (r-r.sub.z) with 45.degree.<.theta..sub.l 135.degree. . All following considerations relate to B.sub.LG,1 (r-r.sub.z); for B.sub.LG,2 (r-r.sub.z) the results are similar.

With r.sub.z .multidot.n.sub.l =x.sub.z cos .theta..sub.l +y.sub.z sin .theta..sub.l and ##EQU20## the following holds: ##EQU21## Likewise, with r.multidot.n.sub.l =x cos .theta..sub.l+y sin .theta..sub.l and ##EQU22## one can write: ##EQU23## It thus follows for B.sub.LG,1 (.rho..sub.x,.rho..sub.y): ##EQU24## is the one-dimensional Fourier transform of the f.sub.lk to the grid ##EQU25## In the frequency space, ##EQU26## is multiplied by the Fourier transform of the convolution Kernel, ##EQU27## and the phase factor ##EQU28## which controls the position of the reconstruction center. There then follows an FFT back into the spatial domain; this results in the convoluted projections f'.sub.c (p.sub.q, .theta..sub.l), centered around -r.sub.z : ##EQU29## In the spatial domain, the f'.sub.c (p.sub.q, .theta..sub.l) are limited to the desired image field region v.sub.l =2D.sub.B cos .theta..sub.l, i.e. to ##EQU30## values, through multiplication by R.sub.l (p.sub.q), and are newly transformed in the frequency space to the grid ##EQU31## with a chirp-z transformation ##EQU32## The spectrum B.sub.LG,1 (.rho..sub.x, .rho..sub.y) is defined in the .rho..sub.x direction at discrete grid points s.DELTA..rho..sub.x.sup.LG, with ##EQU33## B.sub.LG,1 (.sub.x,.rho..sub.y) belongs to an image limited in the x direction with an ideal rectangular window width 2D.sub.B. The grid values s.DELTA..rho..sub.x.sup.LG can thus be used directly as starting points for the FFT in the x direction back into the spatial domain, without requiring an additional interpolation in the frequency space. As in the gridding method, a region of side length D.sub.R =2D.sub.B is thereby reconstructed. If the image D.sub.B .multidot.D.sub.B, having 512.multidot.512 pixels is thereby represented, 1024 frequency values are necessary in the .rho..sub.x direction. If more frequency support points are present in the .rho..sub.x direction after the chirp-z transformation of the projections, because the maximum frequency contained in the convoluted projections is larger than the frequency that can be represented with the image matrix, the additional frequency points must be rearranged as aliasing components after the gridding in the .rho..sub.y direction (see below).

In the .rho..sub.y direction, the grid values are likewise equidistant; the grid depending on s.DELTA..rho..sub.x.sup.LG is ##EQU34## As in the original linogram method, for each s, a chirp-z transformation to the desired pixel grid .DELTA.y can now follow in the .rho..sub.y direction.

In order to avoid the costly chirp-z transformation with the grid .DELTA..rho..sub.y.sup.LG (s), which depends on s.DELTA..rho..sub.x.sup.LG, it is desirable to carry out a one-dimensional interpolation in each column s of the frequency matrix in such a way that a frequency grid .DELTA..rho..sub.y that is equal for all s results, i.e. a Cartesian matrix. In order for no image deterioration is to appear through this interpolation, it must be implemented as a one-dimensional "gridding." The basic idea for the choice of the interpolation function T(.rho..sub.y) is to limit the image B.sub.LG,1 (x,y) to the desired width D.sub.R conceptually through multiplication by a window function T(y) in the y direction, so that no coverage errors occur during the periodic repetition of the image in the spatial domain caused by the spectrum sampling in the Cartesian coordinates. The limiting is not carried out in the spatial domain, but rather in the frequency space, as one-dimensional convolution in the .rho..sub.y direction of the image spectrum defined at discrete points ##EQU35## with the Fourier-transform T(.rho..sub.y) of T(y). The address calculations may thereby be simplified due to the equidistant starting points s.DELTA..rho..sub.y.sup.LG (s). The resulting spectrum, continuous in the .rho..sub.y direction, is then sampled at the grid points j.DELTA..rho..sub.y, with .DELTA..rho..sub.y =1/D.sub.R. In order to obtain the smallest possible extension of T(.rho.) in the frequency space, and thereby as short a convolution as possible, the side length of the reconstruction field D.sub.R is chosen twice as large as the side length of the desired image field D.sub.B. T(y) can then fall gently to a sufficiently small value in the area D.sub.B /2.ltoreq..vertline.y.vertline..ltoreq.D.sub.R -D.sub.B /2.

In the .rho..sub.x direction, no gridding is necessary, since here the Cartesian grid is already present. The image resulting from this spectrum through Cartesian sampling and two-dimensional Fourier back-transformation is identical with the convolution back-projection image in the desired image area up to aliasing errors in the spatial domain, which are caused by the finite extension of T in the frequency space. Through suitable choice of T, these aliasing errors can be made negligibly small.

Through the convolution with T(.rho..sub.y) in the frequency space, the image is multiplied by T(y) in the spatial domain, and thus has a "bump". This bump can easily be corrected by dividing the reconstructed image B.sub.LG,1 (x,y) in the region .vertline.y.vertline..ltoreq.D.sub.B /2 by T(y). After the gridding and a possible rearrangement of the aliasing components, there are 1024 Cartesian values both in the .rho..sub.x direction and the .rho..sub.y direction, so that one obtains the partial image B.sub.LG,1 (r-r.sub.z) through two-dimensional back-transformation of 1024.multidot.1024 frequency points. The inner 512 .multidot.512 pixels are the image that must be corrected in the y direction through multiplication by ##EQU36## In the same way, but with x and y exchanged, the partial image B.sub.LG,2 (r-r.sub.z) is obtained from ##EQU37## Since this partial image must be corrected after the gridding and the back-transformation in the x direction, the spectra of B.sub.LG,1 (r-r.sub.z) and B.sub.LG,2 (r-r.sub.z) cannot simply be added and back-transformed together. Rather, each spectrum must be separately transformed into the spatial domain. The final image results from the sum of the two corrected partial images.

A summary of the invention processing step is as follows:

In the following, the necessary processing steps are shown for a third-generation CT apparatus with 1024 measurement values per projection and .about.1200 fan projections in the region [0,2.pi.] (e.g. Siemens SOMATOM AR/HP).

1.) Rebinning to parallel data f.sub.lk :

N=1024 channels per projection, N.sub.p .apprxeq.600 projections with .DELTA.tan .theta.=const in the angle area [-45.degree., 45.degree.] or .DELTA.cot .theta.=const in the angle area [-45.degree., 135.degree.].

2.) Zero padding of the f.sub.lk and one-dimensional FFT of the f.sub.lk with 2048 values. Through the zero padding, the grid ##EQU38## results in the frequency space.

3.) Multiplication of the f.sub.l (m.DELTA..rho.) with the Fourier transform of the convolution Kernel L.sub.c (m.DELTA..rho.) and with the phase factor exp (-2.pi.im.DELTA..rho.(x.sub.z cos .theta.+y.sub.z sin .theta..sub.l)).

4.) FFT of length 2048 back into the spatial domain; limitation of each convoluted, centered projection f'.sub.c (p.sub.q, .theta..sub.l) to ##EQU39## values for -45.degree..ltoreq..theta..sub.l .ltoreq.45.degree., or to ##EQU40## values for 45.degree.<.theta..sub.l .ltoreq.135.degree..

5.) One -dimensional chirp-z transformation of the bounded function f'.sub.c (p.sub.q,.theta..sub.l) to a maximum of 1024 frequency points in the grid ##EQU41## The number of points depends on the size of the represented image segment and on the maximum frequency of the reconstruction.

6.) Multiplication of the F.sub.l (s.DELTA..rho..sub.F.sup.l) by the Fourier transform of the "smoothing Kernel" L.sub.s (s.DELTA..rho..sub.F.sup.l).

7.) Calculation of two sub-images for the angle areas [-45.degree., 45.degree.] and [45.degree., 135.degree.]: Sub-image 1: for each frequency value s.DELTA..rho..sub.x in the .rho..sub.x direction, a gridding is carried out in the .rho..sub.y direction, on the basis of N.sub.P frequency support points in the s-dependent grid .DELTA..rho..sub.y.sup.LG (s), to 1024 values in the grid ##EQU42## Since the matrix is complexly conjugate, only the positive s must be taken into account; at most there are thus .about.1024.multidot.600 starting points. After the rearrangement of the aliasing components, a Fourier transformation of the 1024.multidot.1024 frequency values into the spatial domain follows; the inner 512.multidot.512 pixels are the image. Sub-image 2: proceed as for sub-image 1, with x and y exchanged.

8.) Correction and addition of the two sub-images.

Although modifications and changes may be suggested by those skilled in the art, it is the intention of the inventors to embody within the patent arranged hereon all changes and modifications as reasonably and properly come with the scope of their contribution to the art.

Claims

1. A method for obtaining an image of a subject in a computed tomography apparatus comprising the steps of:

irradiating a subject with x-rays from a plurality of projections p.sub.k at respective projection angles .theta..sub.l ;

detecting said x-rays, attenuated by said subject, in a radiation detector having N detector channels to obtain data f(p.sub.k, .theta..sub.l);

supplying said data f(p.sub.k, .theta..sub.l) from the detector to a computer in a parallel beam geometry so that .DELTA.tan .theta.=const in an angle region [-45.degree., 45.degree.] or .DELTA.cot .theta.=const in an angle region [45.degree., 135.degree.], and p.sub.k =(k-0.5+a.sub.m)a (a is a sampling grid and a.sub.m is an alignment factor);

reconstructing a sectional image of said subject by the following steps:

one-dimensionally Fast Fourier transforming f(p.sub.k, .theta..sub.l) relative to p.sub.k, with N.sub.FFT .gtoreq.2N values to obtain f.sub.l (m.DELTA..rho.);

multiplying f.sub.l (m.DELTA..rho.) by the Fourier transform of a convolution Kernel L.sub.c (m.DELTA..rho.) and by a phase factor exp (-2.pi.im.DELTA..rho.(x.sub.z cos .theta..sub.l +y.sub.z sin .theta..sub.l)), which defines the position of the center of the sectional image to obtain frequency values on a polar grid;

conducting a Fast Fourier transformation of length N.sub.FFT of said frequency values on a polar grid back into the spatial domain to obtain a plurality of convoluted, centered projections f'.sub.c (p.sub.q,.theta..sub.l);

limiting each convoluted, centered projection f'.sub.c (p.sub.q,.theta..sub.l)to ##EQU43## values for -45.degree..ltoreq..theta..sub.i .ltoreq.45.degree. or to ##EQU44## values for 45.degree.<.theta..sub.i .ltoreq.135.degree. to obtain a plurality of bounded projections f'.sub.c (p.sub.q,.theta..sub.l);

one-dimensionally chirp-z transforming each bounded projection f'.sub.c (p.sub.q, .theta..sub.l) to a grid ##EQU45## for -45.degree..ltoreq..theta..sub.l .ltoreq.45.degree. or ##EQU46## for 45.degree.<.theta..sub.l .ltoreq.135.degree., whereby the number of points is dependent on the size of the represented image segment and on the maximum frequency of reconstruction to obtain projections F.sub.l (s.DELTA..rho..sub.F.sup.l) by the Fourier transform of a smoothing Kernel L.sub.s (s.DELTA..rho..sub.F.sup.l);

calculating first and second sub-images respectively for the angle regions [-45.degree., 45.degree.] and [45.degree., 135.degree.];

calculating said first sub-image by, for each frequency values s.DELTA..rho..sub.x in the p.sub.x direction, interpolating said values with T(.rho..sub.y) in the .rho..sub.y direction, using the gridding method on the basis of N.sub.P frequency support points in an s-dependent grid .DELTA..rho..sub.y.sup.LG (s), to 2N.sub.Pix values in a grid ##EQU47## (N.sub.Pix) is a number of pixels in one column of the sectional image) rearranging aliasing components, and Fourier transforming the 2N.sub.Pix .multidot.2N.sub.Pix frequency values in the spatial domain and correcting inner N.sub.Pix .multidot.N.sub.Pix pixels in the y direction by multiplication with 1/T(y);

calculating said second sub-image by, for each frequency value s.DELTA..rho..sub.y in the .rho..sub.y direction, interpolating said values with T(.rho..sub.x) in the .rho..sub.x direction, using the gridding method on the basis of N.sub.P frequency support points in an s-dependent grid .DELTA..rho..sub.x.sup.LG (s), to 2N.sub.Pix values in a grid ##EQU48## (N.sub.Pix is a number of pixels in one column of the sectional image) rearranging aliasing components, and Fourier transforming the 2N.sub.Pix .multidot.2N.sub.Pix frequency values in the spatial domain and correcting inner N.sub.Pix .multidot.N.sub.Pix pixels in the x direction by multiplication with 1/T(x); and

adding said first and second sub-images.