Iterative Optimizing Quantization Method for Reconstructing Three-Dimensional Images From a Limited Number of Views

1. Field of the Invention

The present invention relates to tomography and more particularly to reconstructing three-dimensional images from a limited number of two-dimensional projections.

2. Description of Related Art

Three-dimensional objects that lie buried away from direct viewing inside a body can often be scanned just outside the body from a two-dimensional plane. Such viewpoints are usually limited in their number by practical considerations. But in general, the more different two-dimensional views that are available, the better the image reconstruction can be of the buried three-dimensional object.

Various methods have been devised to solve the inverse problem for reconstruction with a limited number of views. For example see, R. Gordon, R. Bender and G. T. Herman, "Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography," J. Theor. Biol. 29, 471-481, 1970; P. Gilbert, "Iterative methods for the three-dimensional reconstruction of an object from projections," J. Theor. Biol., vol. 72, pp. 105-117, 1972; A. V. Lakshminarayanan and A. Lent, "Methods of least squares and SIRT in reconstruction," J. Theor. Biol., vol. 76, pp. 267-295, 1979; G. T. Herman, A. Lent, and S. W. Rowland, "ART: Mathematics and applications. A report on the mathematical foundations and on the applicability to real data of the algebraic reconstruction techniques," J. Theor. Biol. vol. 42, pp. 1-32, 1973; R. Gordon and G. T. Herman, "Reconstruction of pictures from their projections," Comm. A.C.M. vol. 14, pp. 759-768, 1971; H. P. Weiss and J. A. Stein, "The effect of interpolation in CT image reconstruction," Proc. IEEE Comput. Soc. Conf. Pattern Recognition Image Processing Chicago, Ill., p. 193, 1978; G. T. Herman, S. W. Rowland, and M. M. Yau, "A comparative study of the use of linear and modified cubic spline interpolation for image reconstruction," IEEE Trans. Nucl. Sci. NS-26, 2879-2894, 1979; and, R. S. Fager, K. V. Peddanarappagari, and G. N. Kumar. "Pixel-based reconstruction (PBR) promising simultaneous techniques for CT reconstructions," IEEE Transactions on Medical Imaging, vol. 12, no. 1, March 1993.

One of the better developed techniques is the algebraic reconstruction technique (ART), introduced by Gordon, et al., which is an iterative procedure for recovering objects from their projections. ART takes a ray at a time and corrects the pixels dynamically as the errors are calculated. But, no averaging is considered. A modification of this method is the simultaneous iterative reconstruction technique (SIRT), proposed by Gilbert in 1972, which has the advantage of reduced noise levels due to averaging but the procedure does not always converge. This divergence was demonstrated by Lakshminarayanan and Lent who proposed a new method called Generalized SIRT (GSIRT), which uses "pseudo-projections" of a fictitious object. In 1971, Gordon and Herman exploited a Monte Carlo method to optimize the image reconstruction and which quantized a digital picture into its picture elements. They demonstrated how to optimize a reconstructed image without disturbing the original projections by adding or subtracting a unit, e.g., quantized, element in each pixel, and also showed the advantage of averaging the signal values for each pixel to enhance the image.

An object of the present invention is to provide a method for reconstructing three-dimensional images from a limited number of two-dimensional projections.

A further object of the present invention is to provide a method for reconstructing three-dimensional images from two-dimensional projections that reduces the artifacts and image distortion that usually accompany input data obtained from a limited number of such two-dimensional projections.

A still further object of the present invention is to provide a method for reconstructing three-dimensional images from two-dimensional projections that simplifies experimental set-ups and supports non-destructive imaging of microstructures.

Briefly, a three-dimensional image reconstruction method embodiment of the present invention comprises treating the object of interest as a group of elements with a size that is determined by the resolution of the projection data, e.g., as determined by the size of each pixel. One of the projections is used as a reference projection. A fictitious object is arbitrarily defined that is constrained by such reference projection. The method modifies the known structure of the fictitious object by comparing and optimizing its four projections to those of the unknown structure of the real object and continues to iterate until the optimization is limited by the residual sum of background noise. The method is composed of several sub-processes that acquire four projections from the real data and the fictitious object, generate an arbitrary distribution to define the fictitious object, optimize the four projections, generate a new distribution for the fictitious object, and enhance the reconstructed image. The sub-process for the acquisition of the four projections from the input real data is simply the function of acquiring the four projections from the data of the transmitted intensity. The transmitted intensity represents the density distribution, that is, the distribution of absorption coefficients through the object.

An advantage of the present invention is that a power window safety device is provided that is inexpensive and reliable.

Another advantage of the present invention is that a sensor is provided that reverses an automatic closing device when solid obstacles are detected in the closure area.

A further advantage of the present invention is that a timing generator is provided that can detect objects that come into close proximity with glass windows.

FIG. 1 is a flowchart of a method embodiment of the present invention;

FIG. 2 is a perspective diagram representing the relation between the intensity I and the mass profile P of a three-dimensional object with a two-dimensional projection as related to the method of FIG. 1;

FIG. 3 is a diagram representing the determination of a pixel position with two projection coordinates as related to the method of FIG. 1;

FIG. 4 represents data used by the method of FIG. 1 in the reconstruction of a 13.times.13 binary object; (a) introduces numerically constructed two-dimensional unknown real and a known arbitrary object, and demonstrates the progress in reconstruction after each loop of the optimization procedure in the first iteration; and

FIG. 5 diagrams a two-dimensional bound region of a primary fictitious object. A reference one-dimensional projection is converged into a region with boundaries that are defined by the intersection of the available views of the real objects. A three-dimensional bound region is similarly reconstructed from two-dimensional projections.

FIG. 1 illustrates an iterative optimizing quantization technique (IOQT) which is a method of the present invention for dealing with a limited number of views, and referred to herein by the general reference numeral 10. The IOQT method 10 includes the optimization of pseudo-projections, as in GSIRT, and the quantization of projections, as in the Monte Carlo method. The IOQT method 10 treats a target object as a group of elements whose size is determined by the resolution of the projection data. The resolution of the projection data is determined, in turn, by the size of the pixels. By using one of the projections as a reference projection, a fictitious object is arbitrarily defined that is constrained by this reference projection. IOQT modifies the known structure of the fictitious object by comparing and optimizing its four projections to those of the unknown structure of the real object and continues to iterate until the optimization is limited by the residual sum of background noise.

The IOQT method 10 comprises a process 12 to acquire a finite number, e.g., four projections from the real data obtained from observations of the target object, a process 14 to generate an arbitrary distribution to define the fictitious object, a process 16 to optimize the four projections, a process 18 to generate a new distribution for the fictitious object, and a process 20 to enhance the reconstructed image. The process 12 is simply a function of acquiring the four projections from the data of the transmitted intensity. The transmitted intensity represents the density distribution, i.e., the distribution of absorption coefficients through the object.

As shown in FIG. 2, the relation between a transmitted intensity, I, and an initial intensity, I.sub.0, as a function of the direction of the projection, .xi..sub..theta., and a two-dimensional density distribution function of the object, .function.(x.sub..theta.,.xi.), can be represented mathematically by,

where .xi..sub..theta. is the axis of distribution for the projection in the .theta.-direction. This can also be written in the form of a mass distribution function in the .theta.-direction, p(.xi..sub..theta.), as ##EQU1## This mass distribution function is used to generate an initial fictitious object. Since the total mass of each slice should be conserved in each direction especially for an experiment, the mass distribution function is normalized according to, ##EQU2## where .sub.rPn (.xi..sub..theta..spsp.1) is the normalized projection data in the direction of .theta..sub.1 and for the r-th slice. When using integer arithmetic, the signal intensifies are amplified by multiplying with a uniform constant to minimize the residual error, e.g.,

where m.sub.r is the total mass of the r-th slice and .function..sub.m is a magnification factor.

To generate an arbitrary distribution for the fictitious object, one of the four projections in the direction defined as .theta..sub.1 is chosen as a reference projection. The remaining pseudo-projections of the fictitious object are then obtained and their relationship to the projections from the real object are

where the subscripts, o and f, represent the real object and the fictitious object respectively and the subscript r denoting the r-th slice is eliminated for simplicity.

The optimization process 16 reorganizes the group of fictitious elements by optimizing the projections of the fictitious object to approach that of the real object. The optimization process 16 is simply to add or subtract a unit element in each pixel which is a quantized intensity equal to 1% of the value of the algebraic multiplication factor .function..sub.m. This procedure includes a subroutine 22 to decide the pixels on which a quantum is added or subtracted, a subroutine 24 to select the pixels by sustaining the chosen constraint, a subroutine 26 to redistribute the fictitious object, and a subroutine 28 to project the fictitious object onto the original view directions. The entire optimization process 16 typically requires looping through these four subroutines several times until reaching a user specified termination criterion.

Some functions for the mathematical interpretation of the first subroutine of the optimization process 16 need defining. In terms of the rotational matrix R(.theta.), the rotated projection is defined as an operator like,

where P.sub.0 is the unrelated projection, that is, ##EQU3## which represents the j-th row projection in the original coordinates. The data is projected in an arbitrary direction .theta. according to,

where .function.(i,j) is the two-dimensional density distribution function of the fictitious object in the original (i,j) coordinates. The difference is expressed between the real and the pseudo projections in the direction of .theta..sub.1 as,

Referring to FIG. 3, since the position of an individual pixel is determined with two projection coordinates, the projection data can be investigated from the four views to assign a probability value to each pixel in the two-dimensional fictitious object. The value of the probability is determined according to the decision process. If .DELTA.d.sup.j' (.xi..sub..theta..spsp.2).ltoreq.0 and .DELTA.d.sup.j" (.xi..sub..theta..spsp.2).ltoreq.0, then the (i.sub.0,j.sub.0)-pixel, which is determined by projection coordinates, (j,j'), does have a high probability for a unit element .delta.(i.sub.0,j.sub.0) to be subtracted. The relation of (i.sub.0,j.sub.0) and (j,j') is deduced as, ##EQU4## from j'=i.sub.0 cos .theta..sub.1 +j.sub.0 sin .theta..sub.1 and j"=i.sub.0 cos .theta..sub.2 +j.sub.0 sin .theta..sub.2. In contrast, if .DELTA.d.sup.j' (.xi..sub..theta..spsp.1).gtoreq.0 and .DELTA.d.sup.j" (.xi..sub..theta..spsp.2).gtoreq.0, then the (i.sub.0, j.sub.0)-pixel does have a high probability for .delta.(i.sub.0,j.sub.0) to be added to the function .function.(i,j). Otherwise the pixels are defined with low probability, which means there is neither an addition nor a subtraction process. After determining the probabilities of all pixels in the two-dimensional object with all combinations of (j,j'), the pixels are selected on which a unit element will be added or subtracted.

The subroutine 24 selects the points for each row which should be added or removed to satisfy the constraints. Due to the constraint d.sub..function. (.xi..sub..theta..spsp.1), the number of .delta..sub.+ (i.sub.0,j=.xi..sub..theta..spsp.1) is equal to that of .delta..sub.- (i.sub.0, j=.xi..sub..theta..spsp.1) for the j-th row. Therefore, an equal number of .delta.+ and .delta..sub.- will be simultaneously added to or subtracted from .function.(i,j) during the redistribution in subroutine 26 of the fictitious object plane in the third subroutine. The subroutine 28 projects the redistributed fictitious object onto the original view directions to generate the mass profiles for use in the next optimization loop that will subsequently use different projection combinations. Addition and subtraction of unit elements are also done for the other combination of projections. The number of optimization procedures is generalized with combination, .sub.n C.sub.m, where n is the total number of projections minus one for the reference projection and m is a smaller number than n. Here, the number of procedures is three.

The image reconstruction is enhanced in the process 20 before a next iteration or prior to the acquisition of the final result. The main purpose of this process is to accelerate the convergence of the reconstruction. A variety of routines can be applied, e.g., a filtering function, a removing of threshold values to reduce unwanted noise, or a smoothing function like convolution is applied. For the application of the threshold procedure, a threshold value is chosen. Since the mass quanta were enhanced with the magnification factor .function..sub.m in order to amplify the signal intensity, each pixel of the image should have at least the value of .function..sub.m. The pixels which have a lower value than that of .function..sub.m are not the ones which contain mass quanta. However, since image reconstruction is never perfect using a limited number of views, a more representative estimate of the threshold value that has been chosen to use in actual simulations is 50% of .function..sub.m.

Another important function of the enhancement process for simulations is to smooth the estimated image. Since the object is defined to be a group of discrete elements and the point spread function of the zone plate is a Gaussian profile, the reconstruction is smoothed by taking the local average value of neighboring pixels. A convolution routine is used with a two-dimensional Gaussian function. For a next iteration, each row of the reconstruction is equalized to sustain the imposed constraints with the equation, ##EQU5## where I.sup.j (i,j) and F.sup.j (i,j) are the extracted j-th equalized and non-equalized reconstructions respectively.

To completely solve a linear equation in the n.times.n unknowns, n.times.n linear equations are needed, n projection data in n columns. But this is impossible in real experiments. With m projections, m<<n, there will be more than one picture with exactly the same projections. However, the process 16 can select points with high probability to be added or removed by comparing the difference between the real and pseudo projections but not those with low probability. After enhancing the reconstruction after each iteration, the probabilities attached to each element will be affected, e.g., some are increased in probability. At this point, those elements will be selected to be added or removed within the fictitious object plane. The process 16 selects the most probable elements by changing their probabilities without disturbing the reference projection. The optimization process 16 continues for the remaining pseudo-projections up until the reconstruction quality degrades. The IOQT method 10 is then terminated. The standard deviation is applied, e.g., the root mean square, of the projection data, given in equation (13), for the m-th loop to terminate each iteration and also for q-th iteration to terminate. For example, ##EQU6## See, W. S. Haddad, I. McNulty, J. E. Trebes, E. H. Anderson, R. A. Levesque, and L. Yang, "Ultrahigh-resolution x-ray tomography," Science, vol. 266, pp. 1213-1215, 1994. According to such termination criterion, if the sequential deviation, .sub.m.sup.q .DELTA., is smaller than the subsequent deviation, .sub.m+q.sup.1 .DELTA., in the q-th iteration, then the loop of the optimization is stopped to generate the q-th estimated image. If the last deviation in q-th iteration, .sup.q .DELTA., is smaller than the subsequent deviation, .sup.q+1 .DELTA., the execution is terminated and the final reconstruction obtained.

The IOQT method 10 was used for the reconstruction of three simulated objects and one experimental test pattern. The results demonstrated the IOQT method 10 is capable of good three-dimensional reconstruction with only a limited number of views. In a first simulation, a 13.times.13 binary image was composed of point objects. FIG. 4 represents how an arbitrarily distributed fictitious object can be quickly massaged to approach a real object's distribution during only three loops of the three combinations of the optimizing procedure over a complete iteration. For convenience, the projection angles used in the test case were zero, forty-five, ninety, and 135 degrees. The projections of the real objects are considered to be real experimental projections, and the projections of the fictitious object are considered to be pseudo-projections.

Two other simulations provided a more rigorous test of the performance of the IOQT method 10. More realistic examples were investigated that had arbitrary projections, gray scale objects, and noise problems. For both the test cases, a 50.times.50 binary image was generated as the real object. In these examples, the finite dimensions of a pixel were taken into account. A value of 100 was used for the magnification factor, .function..sub.m ; a value of fifty was used for the threshold value; and a value of one was used for the standard deviation of the smoothing Gaussian function.

A good agreement was observed between the reconstruction and the real object. The internal structure was revealed, but the choice of projection angles can create a discrepancy in the reconstruction due to the difficult geometries involved. With zero, forty-five, ninety, and 135 degree projections, the reconstruction was good and showed the internal structure. The test object had a symmetry which was projected in evenly separated directions. With zero, fifteen, ninety, and 105 degree projections, artifacts appeared due to the selection of projection angles for an object which was composed of a large number of pixels. The artifacts could be removed by modifying the signal enhancement process.

An improvement was made possible through an improved enhancement procedure that allowed for variations in the threshold value during the image reconstruction for each iteration. In this case, a threshold value of 20% was initially applied for the image reconstruction in the first iteration and then increased by 10% for each subsequent iteration due to the improvement of image perfection. The only difference between the computer simulations was the threshold value that was applied to reduce the residual background. The criterion of equation (13) was used to determine when the IOQT method 10 should be terminated. The simulations demonstrated how the method might be used to select optimum viewing angles. In addition, the optimum signal enhancement technique also appears to depend on the object.

Noise, in the form of random errors like the photon fluctuations produced by unpredictable and unknown sources, is important to consider in simulations. Reference projections that have been contaminated by noise can lead to reconstructions that do not faithfully reproduce the unknown real object. To test the ability of the IOQT method 10 to overcome this problem, a Gaussian random noise distribution was included in each pixel of the simulated projections. The mean value is the noise-free value of each pixel in the simulated projections, and the standard deviation is 10% of the mean value. Two pixels from different projections, representing the position of a point on a two-dimensional object, cannot have the same relative, e.g., corresponding, amount of random noise. As a consequence and after optimization, the noise simply becomes part of the residual background and does not affect the quality of the final image reconstruction. It is assumed, therefore, that the IOQT method 10 is not noise sensitive.

The IOQT method 10 was also tested with a gray scale object in the reconstruction of a real object made of a variety of materials. A gray-scaled object of 50.times.50 pixels was built of several smaller square blocks each of 5.times.5 pixels. The individual pixel intensity within a square block ranged from one to seven. All of the pixels in any given square block had the same value. The intensity of adjacent square blocks differed by one unit. The image reconstructions corresponded to the same set of projection angles mentioned earlier. The differences in the shape and size of the distribution of pixel intensities within the reconstructed images of the square blocks reflected the differences in the value of the initial mass, e.g., absorption coefficient, intensities that were used. The higher the initial mass intensity used, the closer the shape and size of the reconstructed image to the original square block. The simulated real projection of the object also included Gaussian random noise.

A real world three-dimensional reconstruction with the IOQT method 10 used the experimental data collected by Haddad, et al., with a scanning x-ray transmission microscope at Brookhaven National Laboratory. One-dimensional projection data was extracted from each projection with angles of zero, fifteen, and .+-.55 degrees for one of a two-dimensional object plane or slice. In a series of reconstructions, each iteration improved on removing artifacts and reducing the residual background until the real test pattern was approximated.

A critical part of any optimization method is the selection of a convergence criterion. The reconstruction cannot be uniquely determined by a finite number of projections. The criterion for the IOQT method 10 has been found to be consistent with the terminating criterion defined by equation (13). The reconstructed images have all been found to converge into a region with boundaries defined by the intersection of the four views of the real object. Different combinations of primary distributions for the fictitious object and reference projections were used as constraints. The differences between the reconstructed images and actual were indistinguishable. If a reconstructed image did not fall within such boundary, then the reconstructed image did not represent the information obtained by all the limited number of views of the real object. This is exactly the condition that would be required for the standard deviation to be at a minimum.

An example for the highly probable region for solution convergence in the case of four views using the IOQT method 10 is illustrated in FIG. 5. In each iteration, the reconstruction with the IOQT method 10 is always convergent. The optimization procedure 16 minimizes the difference, .DELTA.d.sup.j =d.sub.0.sup.j -d.sub.f.sup.j. The estimated image is bounded by a plausible region. However, the sequential ratio for the q-th iteration, ##EQU7## is not always less than one. Since the signal enhancement procedure redistributes the intensity profile of the reconstruction as described by equation (12) in order to facilitate the next iteration, the ratio of standard deviations is greater than one. Then IOQT method 10 is terminated. Averaging the reconstructed images obtained from using different combinations of the primary distribution of the fictitious object and reference projections, terminated using the criterion in equation (13), would lead to even a better image by reducing the overall number of artifacts.

The present invention thus includes a method of reconstructing a three-dimensional image in a volume from a finite number of two-dimensional projections at various limits of the volume. The method comprises collecting a plurality of two-dimensional projections of a three-dimensional object-of-interest in a volume from various points of observation outside the volume: (1) Defining the three-dimensional object-of-interest as a group of elements; (2) selecting any one of the plurality of two-dimensional projections to serve as a reference projection; (3) arbitrarily defining a fictitious three-dimensional strawman in the volume that is constrained by the reference two-dimensional projection; (4) modifying the three-dimensional strawman one at a time to conform to each of the other two-dimensional projections; (5) optimizing by continued iterations through the step of modifying and comparing until further optimization of the three-dimensional strawman is limited by a residual sum of background noise; and (6) outputting a representation of the optimized strawman as an estimate of the three-dimensional character of the object-of-interest in the volume.

The defining of the three-dimensional object-of-interest in the volume as the grouping of elements can be such that each element has a size that depends on the resolution of a particular two-dimensional projection data collected, as is ordinarily determined by the size of corresponding image pixels.

The step of collecting can include the selection of a finite number of projections from an observation data input that is a function of acquiring such projections from transmitted intensity data, wherein the transmitted intensity represents a density distribution of absorption coefficients through the three-dimensional object-of-interest.

Thus another method of reconstructing a three-dimensional image in a volume from a finite number of two-dimensional projections at various limits of the volume comprises acquiring a fixed number of two-dimensional projections from a fictitious object: (1) Generating an arbitrary distribution to define the fictitious object; (2) optimizing the fixed number of projections; (3) generating a new distribution for the fictitious object; and (4) then enhancing the reconstructed image.

Although particular embodiments of the present invention have been described and illustrated, such is not intended to limit the modifications and changes will no doubt become apparent to those skilled in the art, and it is intended that the invention only be limited by the scope of the appended claims.