Enhanced local tomography

Abstract

Local tomography is enhanced to determine the location and value of a discontinuity between a first internal density of an object and a second density of a region within the object. A beam of radiation is directed in a predetermined pattern through the region of the object containing the discontinuity. Relative attenuation data of the beam is determined within the predetermined pattern having a first data component that includes attenuation data through the region. In a first method for evaluating the value of the discontinuity, the relative attenuation data is inputted to a local tomography function .function..sub..LAMBDA. to define the location S of the density discontinuity. The asymptotic behavior of .function..sub..LAMBDA. is determined in a neighborhood of S, and the value for the discontinuity is estimated from the asymptotic behavior of .function..sub..LAMBDA.. In a second method for evaluating the value of the discontinuity, a gradient value for a mollified local tomography function .gradient..function..sub..LAMBDA..epsilon. (x.sub.ij) is determined along the discontinuity; and the value of the jump of the density across the discontinuity curve (or surface) S is estimated from the gradient values.

Description

BACKGROUND OF THE INVENTION

This invention relates to image reconstruction from tomographic data, and, more particularly, to the definition of discontinuity location and size using limited tomographic data. This invention was made with government support under Contract No. W-7405-ENG-36 awarded by the U.S. Department of Energy. The government has certain rights in the invention.

Tomography produces the reconstruction of a generalized density function f from a large number of line integrals of f. A practically important objective is the reconstruction, from the x-ray data, of functions providing significant information about the physical object, such as the location of discontinuities within an object being interrogated. Tomography equipment for obtaining the attenuation data is well known. For example, in some instances a single source is collimated and traversed across the object, whereby a single sensor output corresponds to a single source transverse location. Here, the source traverses the object at each angular orientation of the object. In other instances, a single source generates a fan-like pattern of radiation that is detected by an array of sensors, where the object is located between the source and the sensor array. The source is then angularly rotated relative to the object to provide a complete set of tomographic data.

Conventional tomography is a global procedure in that the standard convolution formulas for reconstruction of the density at a single point require the line integral data over all lines within some planar cross-section containing the point. Local tomography has been developed for the reconstruction at a point where attenuation data is needed only along lines close to that point within the same cross-section. Local tomography produces the reconstruction of a related function using the square root of the negative Laplace operator and reproduces the locations of discontinuities within an object. See, e.g., E. Vainberg, "Reconstruction of the Internal Three-Dimensional Structure of Objects Based on Real-Time Integral Projections," 17 Sov. J. Nondestr. Test., pp. 415-423 (1981); A. Faridani et al., "Local Tomography," 52 SIAM J. Appl. Math, No. 2, pp. 459-484 (April 1992).

Local tomography can reduce significantly the amount of data needed for the local reconstruction, with a concomitant reduction in x-ray dose. While the location of a discontinuity is reproduced, however, them is no quantitative value for the magnitude of the discontinuity. In many instances it would be useful to know this value in order to make medical, technological, or geophysical conclusions.

Accordingly, it is an object of the present invention to determine both the location and size of discontinuities from tomographic data.

It is another object of the present invention to determine the location and size of discontinuities from only limited attenuation data that includes a region containing the discontinuity.

Additional objects, advantages and novel features of the invention will be set forth in part in the description which follows, and in part will become apparent to those skilled in the art upon examination of the following or may be learned by practice of the invention. The objects and advantages of the invention may be realized and attained by means of the instrumentalities and combinations particularly pointed out in the appended claims.

SUMMARY OF THE INVENTION

To achieve the foregoing and other objects, and in accordance with the purposes of the present invention, as embodied and broadly described herein, the apparatus of this invention may comprise a method for determining by tomography the location and value of a discontinuity between a first internal density of an object and a second density of a region within the object. A beam of radiation is directed in a predetermined pattern through the region of the object containing the discontinuity. Relative attenuation data of the beam is determined within the predetermined pattern having a first data component that includes attenuation data through the region.

In a first method for evaluating the value of the discontinuity, the relative attenuation data is inputted to a local tomography function .function..sub..LAMBDA. to define the location S of the density discontinuity. The asymptotic behavior of .function..sub..LAMBDA. is determined in a neighborhood of S, and the value for the discontinuity is estimated from the asymptotic behavior of .function..sub..LAMBDA..

In a second method for evaluating the value of the discontinuity, a gradient value for a mollified local tomography function .gradient..function..sub..LAMBDA..epsilon. (x.sub.ij) is determined along the discontinuity; and the value of the jump of the density across the discontinuity curve (or surface) S is estimated from the gradient values.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings, which are incorporated in and form a part of the specification, illustrate the embodiments of the present invention and, together with the description, serve to explain the principles of the invention. In the drawings:

FIG. 1 is a graphical representation of the functions .psi.(q).

FIG. 2 is a graphical representation of the functions .psi.'(q).

FIG. 3 illustrates an exemplary phantom for generating Radon transform data.

FIG. 4 is a horizontal cross-section of the density plot through the phantom shown in FIG. 3.

FIG. 5 is a horizontal cross-section of the local tomographic density discontinuity values obtained according to a first reconstruction process of the present invention.

FIG. 6 is a horizontal cross-section of the local tomegraphic density discontinuity values obtained according to a second reconstruction process of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

In accordance with the present invention, local tomography is enhanced to determine the sizes of density discontinuities within an object, as well as the location of the discontinuity. The density function f(x) is usually reproduced from the line integral data .function.(.theta.,p), where .theta. is an angle on the unit circle and p is a location of a beam line, by the formula: ##EQU1##

A local tomography function .function..sub..LAMBDA. is proposed in, e.g., E. Vainberg, "Reconstruction of the Internal Three-Dimensional Structure of Objects Based on Real-Time Integral Projections," 17 Sov. J. Nondestr. Test., pp. 415-423 (1981), and K.T. Smith et al., "Mathematical Foundations of computed Tomography," Appl. Optics, pp. 3950-3957 (1985): ##EQU2## where .THETA..di-elect cons.S.sup.1, .theta. is a scalar variable 0<.theta.<2.pi., .function..sub.pp is the second derivative with respect to p, and S.sup.1 is the unit circle.

Then, it can be shown that .function..sub..LAMBDA. satisfies ##EQU3## where F and F.sup.-1 denote Fourier transform and its inverse, respectively, and (-.DELTA.).sup.1/2 is the square root of the negative Laplacian.

In practice, the function .function..sub..LAMBDA. is not computed, but its mollification W.sub..epsilon. *.function..sub..LAMBDA.. A mollifier is a kernel such that its convolution with f(x) gives a sufficiently smooth function that converges to f(x) as .epsilon..fwdarw.0. Let W.sub..epsilon. be a sequence of sufficiently smooth mollifiers that satisfy the following properties:

(a) W.sub..epsilon. (x) is a radial function, W.sub..epsilon. (x)=W.sub..epsilon. (.vertline.x.vertline.);

(b) W.sub..epsilon. (x)=0, .vertline.x.vertline.>.epsilon.;W.sub..epsilon. (x)>0,.vertline.x.vertline.<.epsilon.; and

(c) W.sub..epsilon. (x)=.epsilon..sup.-2 W.sub.1 (x/.epsilon.), .intg..sub..vertline.x.vertline.<1 W.sub.1 (x)dx=1.

It can be shown that ##EQU4## where R is the Radon transform. Integrating by parts, one gets ##EQU5## From Eq. (3) it is seen that .function..sub..LAMBDA. is the result of the action of pseudo-differential operator (PDO) with symbol .vertline..xi..vertline. on f. Let S be a discontinuity curve of f and pick x.sub.0 .di-elect cons.S, so that S is smooth in the neighborhood of x.sub.0. The asymptotic formula for .function..sub..LAMBDA. in a neighborhood of x.sub.0 can be shown to be ##EQU6##

Combining Eqs. (4) and (5), and using the properties of W.sub..epsilon., for a sufficiently small .epsilon.>0 so that S can be considered flat inside an .epsilon.-neighborhood of x.sub.0 : n.sub.+ (x.sub.S)=n.sub.+ (x.sub.0)=n.sub.+, .vertline.x.sub.S -x.sub.0 .vertline.<.epsilon., yields ##EQU7## where .eta..sub..epsilon. =.eta.*W.sub..epsilon.. Since w.sub..epsilon. =RW.sub..epsilon., ##EQU8##

Since .eta. is smooth in the neighborhood of x.sub.0, the function .eta..sub..epsilon. =.eta.*W.sub..epsilon. is also smooth and variation of .eta..sub..epsilon. is neglected compared with the variation of .psi. on distances on the order of .epsilon. from x.sub.0. Thus, from Eq. (8) and letting q=h/.epsilon., an approximate formula for .function..sub..LAMBDA..epsilon. in a neighborhood of x.sub.0 is: ##EQU9## where c=.eta..sub..epsilon. (x.sub.S). In the case of the mollifier ##EQU10## the graphs of .psi. and .psi.' are shown in FIGS. 1 and 2, respectively.

Equation (9) explicitly relates the behavior of .function..sub..LAMBDA..epsilon. in a neighborhood of x.sub.S with the unknown quantity D(x.sub.S), which is the value of the jump of f at x.sub.S .di-elect cons.S. Since .function..sub..LAMBDA..epsilon. is computed by Equation (4) from the tomographic data and the function .psi.(q) can be tabulated for a given sequence of mollifiers, an estimation of the value of the density discontinuity jump D(x.sub.S) of the original function f at x.sub.S is provided by Equation (9).

A first process for evaluating D(x.sub.S) begins with the computed local tomography function calculated on a square grid with step size h:x.sub.i,j =(x.sub.i.sup.(1),x.sub.j.sup.(2))=(ih,jh), where i and j are integers. Choose a grid node x.sub.i.sbsb.0.spsb.j.sbsb.0 on S, and assume that h and .epsilon. are sufficiently small that changes of .function..sub..LAMBDA..epsilon. can be neglected in the direction parallel to S, i.e., assume D(x.sub.S)=D(x.sub.i.sbsb.0.spsb.j.sbsb.0) if .vertline.x.sub.S -x.sub.i.sbsb.0 .spsb.j.sbsb.0.vertline.<.epsilon.. Then Eq. (9) can be written as ##EQU11##

Fix n.sub.1, n.sub.2 .di-elect cons.N and consider a (2n.sub.1 +1).times.(2n.sub.2 +1) window around x.sub.i.sbsb.0.spsb.j.sbsb.0. First, estimate n.sub.+ by computing partial derivatives ##EQU12## The estimate is accurate up to a factor (-1), i.e., N.sub.+ =.+-.n.sub.+. Then find D(x.sub.i.sbsb.0.spsb.j.sbsb.0) by solving the minimization problem ##EQU13## Note that the larger values of .function..sub..LAMBDA..epsilon. correspond to the side of S with larger values of f.

Thus, the following process provides estimates of jumps of f from .function..sub..LAMBDA..epsilon. :

1. Estimate vector N.sub.+ from Eq. (13);

Calculation of .vertline.D(x.sub.i.sbsb.0.spsb.j.sbsb.0).vertline. from Eqs. (14) and (15) is most stable if the function .psi. on the interval ##EQU14## differs from a constant as much as possible. For example, the graph of .psi.(q) shown in FIG. 1 suggests that the interval X should be close to the interval between q-coordinates of the local minimum and the local maximum of .psi.(q), i.e., X.apprxeq.[-0.31, 0.31].

A second process for evaluating D(x.sub.S) from the local tomography approach provides values of the jumps from tomographic data in one step. Taking the gradient on both sides of Eq. (12) provides an approximate equality ##EQU15## The modulus of the gradient of .vertline..gradient..function..sub..LAMBDA..epsilon. (x.sub.ij).vertline. can be calculated (Eq. (4)) to get ##EQU16## For each .theta. it is convenient to compute ##EQU17## for all q in the projection of the support of f on the direction .THETA., and then compute ##EQU18## Thus, the following process provides estimated values of jumps of f:

1. Compute .gradient..function..sub..LAMBDA..epsilon. (x.sub.ij) using Eqs. (19), (20), (21);

2. Estimate .vertline.D(x.sub.ij).vertline. using Eqs. (17), (18);

3. The vector .gradient..function..sub..LAMBDA..epsilon. (x.sub.ij) points from the smaller values of f to the larger values of f.

For the mollifier represented by Eq. (11), .psi.'(0) can be computed using Eq. (7): ##EQU19##

FIG. 3 illustrates an exemplary phantom for generating the Radon transform data. The densities are: exterior: 0, ellipse: 1.0, exterior annulus: 0.8, area between the annulus and the ellipse: 0, three small circles: 1.8, center circle: 0.1. The Radon transform was computed for 350 angles equispaced on [0, .pi.] and 601 projections for each angle.

The horizontal cross-section of the density distribution of the phantom is shown in FIG. 4. The function, .function..sub..LAMBDA..epsilon. was computed at the nodes of a square 201.times.201 grid using the mollifier of Eq. (11) with .epsilon.=9.DELTA.p, where .DELTA.p is the discretization step of the p-variable. This value of .epsilon. means that the discrete convolutions were computed using 19 points per integral. The grid step size was h=3.DELTA.p, so that h/.epsilon.=1/3.

A central horizontal cross-section of the density distribution represented by the estimated .vertline.D(x).vertline. from the first process is shown in FIG. 5, with n.sub.1 =n.sub.2 =1, i.e., a 3.times.3 window. It is seen that for the given ratio h/.epsilon.=1/3 and window size n.sub.1 =n.sub.2 =1, the interval X is very close to optimal. The line with peaks is the computed .vertline.D(x).vertline.; the dots represent the actual value of the jumps of the original density function f (see FIG. 4). The estimated jump values are in good agreement with the actual values.

The second process was used on the same phantom structure shown in FIGS. 3 and 4. The grid parameters are the same as for the first process. Again, the line with peaks is the computed .vertline.D(x).vertline.; the dots represent the actual jump values of f. Again, there is good agreement between the actual jump values and the estimated jump values.

The foregoing description of the invention has been presented for purposes of illustration and description and is not intended to be exhaustive or to limit the invention to the precise form disclosed, and obviously many modifications and variations are possible in light of the above teaching. The embodiments were chosen and described in order to best explain the principles of the invention and its practical application to thereby enable others skilled in the art to best utilize the invention in various embodiments and with various modifications as are suited to the particular use contemplated. It is intended that the scope of the invention be defined by the claims appended hereto.

Claims

1. A method for estimating by tomography the location and value of a discontinuity between a first internal density of an object and a second density of a region within said object, comprising the steps of:

directing a beam of radiation in a predetermined pattern through said region of said object containing said discontinuity;

determining relative attenuation data of said beam within said predetermined pattern having a first data component that includes attenuation data through said region;

inputting said relative attenuation data to a local tomography function .function..sub..LAMBDA. to define said location S of said discontinuity;

determining the asymptotic behavior of .function..sub..LAMBDA. in a neighborhood of S; and

estimating said value for said discontinuity from said asymptotic behavior of .function..sub..LAMBDA..

2. A method according to claim 1, wherein said asymptotic behavior of .function..sub..LAMBDA. is determined with the function .psi. computed from a selected sequence of mollifiers, where ##EQU20##

3. A method according to claim 1, wherein said step of determining said asymptotic behavior includes the steps of:

estimating the vector N.sub.+, where ##EQU21## to determine the direction of variation of said local tomography function .function..sub..LAMBDA. about said discontinuity location S;

inputting said vector N.sub.+ to ##EQU22## for estimating said value of said discontinuity.

4. A method according to claim 2, wherein said function .psi. is evaluated over an interval over which .psi. has the greatest variation.

5. A method according to claim 3, wherein said asymptotic behavior of .function..sub..LAMBDA. is determined from ##EQU23## with the function .psi. computed from a selected sequence of mollifiers.

6. A method according to claim 5, wherein said function .psi. is evaluated over an interval over which .psi. has the greatest variation.

7. A method for estimating by tomography the location and value of a density discontinuity between a first internal density of an object and a second density of a region within said object, comprising the steps of:

directing a beam of radiation in a predetermined pattern through said region of said object containing said discontinuity;

determining relative attenuation data of said beam within said predetermined pattern having a first data component that includes attenuation data through said region;

inputting said relative attenuation data to a local tomography function .function..sub..LAMBDA. to define said location S of said discontinuity;

compute a gradient value .gradient..function..sub..LAMBDA..epsilon. (x.sub.ij) along said discontinuity where, ##EQU24## and estimating the value .vertline.D(x.sub.ij).vertline. of said density discontinuity along S by inputting said gradient value to .vertline.D(x.sub.i.sbsb.0.spsb.j.sbsb.0).vertline..apprxeq..vertline..gradient..function..sub..LAMBDA..epsilon. (x.sub.i.sbsb.0.spsb.j.sbsb.0).vertline..pi..epsilon..sup.2 /.psi.'(0).

8. A method according to claim 7, wherein said gradient value indicates the direction of change for said density discontinuity.