Image reconstruction by minimizing second-order moment

ULTRASONIC

IMAGING

IMAGE

12,

71-83

(1990)

RECONSTRUCTION

BY MINIMIZING

Kil-Houm

SECOND-ORDER

Park and Song-Bai

MOMENT

Park

Department of Electrical Engineering Korea Advanced Institute of Science and Technology P.O. Box 150, Chongyangni, Seoul, Korea 131

Statistical methods have received much attention in image reconstruction, especially when the projection data are not sufficient in number, since they usually give more reliable results in such a case than other methods. In this paper, we propose a new statistical method which minimizes the second-order moment of an object that indicates the variance and randomness in a statistical structure. By minimizing the second-order moment, the resulting image reflects the statistical structure imposed by the available information and is biased toward a flat gray structure in the absence of information. The computer simulation results show a good convergence behavior of the proposed method. This method was successfully applied to ultrasound attenuation CT using a sponge phantom. @1990 Academic Press, Inc. Key words

:

A priori information; minimum second-order ment; ultrasound attenuation CT.

moment;

second-order

mo-

I. INTRODUCTION Since the introduction of computed tomography (CT), many image reconstruction methods have been suggested, including statistical methods, which were motivated by the need for obtaining reliable results when the number of projection data are not sufficient. tion, form

An image reconstruction f(z, y), from a finite

m;

= Ill

where

problem is related to estimate a two-dimensional number of linear measured projection data, m;,

i = 1,2,. . . ,h4

hi(s,y)f(r,y)dsdy,

hi(z, y) is the measurement

funcof the

kernel to be specified

depending

on the application.

In many image reconstruction problems, the basic problem is the unavoidable insufficiency of measured projection data, and hence we naturally look for some statistical approaches to the solution. In this case, any finite representation of the object is underdetermined, resulting in an infinity of possible solutions. In fact, the measured data can tell us only that the solution must belong to a specified class of the solutions, but there is no guidance as to what unique choice is to be made within that class. The regularization idea is to recognize the need of employing all of the prior information available, which is of utmost importance in overcoming the problem [I].

71

0161-7346/90 $3.00 Copyright 0 1990 by Academic Press, Inc. All rights of reproduction in any form reserved.

PARK

AND

PARK

In many practical situations, one is able to obtain some information about the object before reconstruction. In addition to the knowledge of the spatial extent of the object, one may also have available a priori information about the range of variation of the physical parameters of interest. The problem is how the a priori information with both of these. The solving inverse problems the a priori information this property into other One statistical

to incorporate the information from the projection data and in the reconstruction procedure to obtain an image consistent statistically motivated approach is a coherent approach for because it allows us to take into account in a natural way in the solution process. In contrast, it is difficult to introduce reconstruction approaches.

approach

is to use the Bayes’ theorem

[2]

which calculates the a posteriori probability P(f/mi) so that any reconstruction f^is corrected by the available measurements. This conditional probability can be maximized if the terms in the numerator of Eq. (2) are once specified; the denominator is a constant and does not affect the optimization. Assuming the measurement equation is known, determination of P(rni/fi is usually straightforward. The difficulty with the Bayes’ theorem is that it depends on the a priori probability density P(T). To find this, there must be a statistical model that gives, for every possible object, an a priori probability of occurrence. Another approach is the maximum entropy reconstruction method [3-161. This method is especially powerful when the projection data are degraded by noise. Enis used to indicate randomness or uncertainty in the tropy Cfl fb, Y) ln[fb y)Wd statistical structure. Projection data in the maximum entropy method reduce randomness and restrict statistical freedom. The principle of the maximum entropy method assures that the description of the situation reflects only the statistical structure imposed by available data. In this paper, we will consider a new approach motivated by non-Bayesian, but still probabilistic arguments. In abstract terms, the procedure begins with the assignment of the second-order moment(SOM) to indicate the randomness and variance in a statistical structure. The concept of the SOM is a familiar one in the probability theory, although its significance in the present problem of image reconstruction is not yet clear. To define the SOM measure that is useful in the context of reconstruction, we must formulate a statistical model for the image processing. Once a given model is adopted, image reconstruction using the SOM proceeds by finding the reconstruction that minimizes the SOM measure and is compatible with all the available information, both projection data and u priori information. Each projection data or piece of a priori information reduces the randomness and restricts the statistical variance. This minimum secondorder moment MSOM) is an intrinsic way of selecting a single image from many images that fit t 6 e available information. The basis of the method is that the resulting image possesses minimum configurational information , so that there must be the evidence in the data for any structure which is seen. If the data are not sufficient MSOM can be used to assure that

to eliminate the resulting

72

all randomness, image becomes

a principle the canonical

of the form.

IMAGE

The

canonical

information

form

reflects

only

RECONSTRUCTION

the statistical

structure

imposed

by the available

and is biased toward a flat gray structure in the absence of information.

Thus, the MSOM description retains all of the randomnessnot removed by the available information, and it is interpreted as the description that is most objective or maximally

noncommittal

II. MINIMUM

with

respect

SECOND-ORDER

to nonavailable

information.

MOMENT

This section describes the procedure of minimizing the SOM. This procedure requires the solution of a constrained optimization problem, to which we will apply the Lagrange multiplier technique[l7]. Under

a statistical

model,

p(x3y)

the probability

P(z, y) of the object

f(?Y) = Js, f(x, y)dxdy

f(X,Y) = F’

where F is the total intensity and D is the cross section the probability distribution M2 is given by M2 =

f 2(x, YP(X,

is given by

(3)

of the object.

The

SOM

of

YWY.

Because this quantity depends not on the total intensity but only on the distribution of gray levels in the image, we will minimize the following quantity :

M; =

(5)

Note that the minimization of Eq. (4) is equivalent to the minimization under the assumption that F is constant, which is valid in the present can be determined uniquely from the projection data. In figure 1, projection be available in the form

Pi,=

J

of Eq. (5) case since F

data Pi,,, of the m-th ray in the j-th view are assumed

dt.f(s,.cosBi-t.sinBi,s,.sinBj+t.cosBj).

In practice, it is known that in eneral there exist satisfy Eq. (6) f or a given set of f Pi,}. We propose select a reliable and safe solution which reflects the the available information and is biased toward a flat information.

73

to

(6)

an infinity of solutions f which the following MSOM method to statistical structure imposed by gray structure in the absence of

PARK

t

AND

PARK

Y

X

% J=l,...,J m=l,..,M

Fig. 1 Parallel

ray geometry.

The MSOM method in image reconstruction intends to minimize Mi in Eq. (5) subject to the constraints given in Eq. (6). By introducing a Lagrange multiplier A to each constraint, we can form the Lagrangian \k as

Q(f,A) =Jl,

f3(s,y)dsdy - k 5 Aim x j=lm=l (7)

[Pi,-

dt.f(s,.cosBj-t.sinBi,s,.sinBi+t.cosBi)], J

where

J and M are the numbers The functional

derivative

f(x, y) = [ 5 2

of views and rays, respectively. of 0 with

2

respect

Ai,,, . xjm(z

to f is set equal to zero, resulting

. cos 0, + y . sin ei)] 1’2,

in

(8)

j=lm=l

where

xi,,,

is the characteristic

xjm(z.cos6i+y.sinBi)

function

=

at the point

S,,,, or

for s, = x.cosBi+y.sinBj otherwise

o1, 3

74

(9)

IMAGE

RECONSTRUCTION

The coefficients {A,,} are determined tions G(A) = 0 obtained by substituting

in principle from the nonlinear Eq. (8) into Eq. (6), or

Ajm

G(A) = Pjm

*Xjm(l.7Z*COS6j

system equa-

+ y*sinBj)]1’2& (10)

= 0.

The image (10).

can be reconstructed

from Eq.

(8) after

A has been obtained

In Eq. (lo), the resulting J x A4 equations in J x it4 unknowns nonlinear, and can be solved, e. g., by an J x M- dimensional nonlinear method, which performs an iterative procedure given by

A(?+‘) =& 1m

where

J(A)

is the Jacobian

matrix

In practice, we do not invert correction term to Ak,

with

-

Eq.

{A .m} are d ewton’s

(11)

J(A”)-‘G(A)

aG(A)/dA.

J(Ak).

J(Ak)Ak+’

from

Instead

= -G(A”).

we solve a linear

system

for a

(12)

After the solution of the linear equations, the new approximation, Ak+‘, are obtained by adding Ak+r to A k. The above Newton’s method is generally expected to give quadratic convergence, provided that a sufficiently accurate starting value is known and J(A) exists.

III. SIMULATION

AND

EXPERIMENT

The MSOM algorithm is especially useful in the situation where the classical linear methods such as convolution backprojection fail to give satisfactory results. When the projection data are sufficient in number to reconstruct an artifact-free image, the use of the MSOM algorithm may not give better results than the linear methods; in addition, its computational cost will be huge. Hence its usefulness may be limited to the situation where the projection data are not sufficient in number.

In the simulations, we consider three cases: 1. 8 views and 65 rays in each view; 2. 16 views and 65 rays; and 3. 24 views and 65 rays.

75

PARK

AND

PARK

1 #8 06 ,4 .2 0. -.2 -.4 -.6 -.8 -I ‘-l-,8-,6-.4-.2

0. ,2 .4

.6

.8

Fig. 2 Test phantom.

1

To compare the results obtained by the proposed algorithm with those from the linear method, we use the convolution backprojection (CB) algorithm. The MSOM algorithm was tested by computer simulation with a test phantom fabricated in our laboratory. As shown in figure 2, the phantom consists of 4 cylindrical regions of the same diameter within a large cylinder. The physical size of each region can be read from figure 2, and the number in each region indicates its attenuation coefficient. All the images were reconstructed on a square pixel of 64 by 64. The image quality is evaluated by the distance defined as

(13) i=lj=l

i=li=l

where f;i is an estimated image at the (ii) pixel, & is the original image, f is the mean of the original image, and Np is the total pixel number of a row or a column. For the test phantom, the distance for the MSOM and the CB are shown in table I. We note that the MSOM algorithm converges to the minimum distances as the iteration proceeds, and the MSOM is superior to the CB especially when the number of view is small. For comparison, the original image is shown in figure 3, and the images reconstructed for the three cases using the MSOM and the CB are shown in figures 4, 5, and 6, where (a) is the reconstructed image by the CB and (b)is the reconstructed image after five iterations by the MSOM. Next, the MSOM algorithm was experimentally tested by ultrasonic attenuation CT. A sponge phantom with five holes, which is widely used in ultrasonic experiments because of its acoustic characteristics are similar to tissue, was scanned by a precisely controlled mechanical scanner. The ultrasonic impulsive wave was generated by a weakly focused 2.25 MHz transducer with a bandwidth of 0.9 MHz. The transmitted wave was received by another transducer having the same characteristics as the transmitter. The test phantom and two transducers were immersed in a water bath and scanned linearly by a stepping motor at 51 linear positions in each of 9, 18, and

76

IMAGE

RECONSTRUCTION

Table I. Distances for the test phantom. Iteratia n

.I = 8, M = 65

J =16, M = 65

1 .1344

J = 24, M = 65

1

1 .1243

1

27 views around the object over 180 degrees. The received signals were digitized a 100 MHz rate by the waveform recoder, and ultrasonic attenuation tomogram reconstructed by the peak value method [18].

Fig.

3 Original

image of the test phantom.

at is

PARK AND PARK

(a)

(b) Fig. 4 Reconstructed reconstructed

images of the test phantom(case by the CB. (b) Image reconstructed

78

1 : J=8, M=65). by the MSOM.

(a) Image

IMAGE

Fig. 5 Reconstructed reconstructed

RECONSTRUCTION

images of the test phantom(case by the CB. (b) Image reconstructed

79

2 : J=16, M=65). by the MSOM.

(a) Image

PARK

AND

PARK

(a)

(b) Fig. 6 Reconstructed reconstructed

images of the test phantom(case by the CB. (b) Image reconstructed

80

3 : J=24, M=65). by the MSOM.

(a) Image

IMAGE

Fig.

RECONSTRUCTION

7 Ultrasonic attenuation tomograms of a sponge(case 1 : J=9, M=Sl). by the MSOM. reconstructed by the CB. (b) I ma g e reconstructed

(a) Image

Figures 7,8, and 9 show the reconstructed images of the sponge phantom by the MSOM and the CB, respectively, with 9, 18, and 27 views. In these figures, (a) shows the reconstructed image by the CB and (b) s h ows the reconstructed images after five iterations by the MSOM. Through simulation and experiment, it is observed that images reconstructed by the MSOM algorithm show substantial reduction of artifacts and better edge conversation compared to those by the CB algorithm. Especially, it is proved that the MSOM reconstructs a better replica of the original image when the projection data are not sufficient in number.

Fig.

8 Ultrasonic attenuation tomograms of a sponge(case reconstructed by the CB. (b) Image reconstructed

81

2 : J=18, M=Sl). by the MSOM.

(a) Image

PARK

AND

PARK

ig. 9 Ultrasonic attenuation tomograms of a sponge(case reconstructed by the CB. (b) Image reconstructed

3 : J=27, M=51). by the MSOM.

(a) Image

IV. CONCLUSION A minimum second-order moment (MSOM) method has been proposed and tested by computer simulation. This method is a new statistical image reconstruction method which minimizes the second-order moment(SOM of an object that indicates the variance or randomness in a statistical structure. alT h e image by the MSOM gorithm reflects the statistical structure imposed by the available information and is biased toward a flat gray structure in the absence of information. Through computer simulation and experiments, it has been shown that the MSOM behaves in a stable manner. In particular, the MSOM is superior to the CB when the projection data are not sufficient in number.

REFERENCES

[II

Sanz, J. L. and Huang, T. S., Discrete and continuous trapolation, IEEE Trane. ASSP 31 1276-1285 (1983).

PI

Richardson, W. H., Bayesian-based Sm. Am. 62, 55-59 (1972).

PI

Frieden, B. R., Restoring with maximum Opt. Sot. Am. 62, 511-518 (1972).

141inWemecke, SPSE PI

iterative

method likelihood

band-limited

signal

of image restoration, and maximum

ex-

J. Opt.

entropy,

J.

S. J., Maximum Entropy Techniques for Digital Image Reconstruction, Conference Proc. R. Shaw, ed., pp. 238243 (1976).

Lent, A., A Convergent Algorithm with a Medical X-ray Application, 249-257 (1976).

S. J. 161Wernecke, Tkans. Computer

and D’Addario, C-26, 351-364

for Maximum Entropy Image Restoration, in SPSE Conference Proc. R. Shaw, ed., pp. L. R., Maximum (1977).

82

entropy

reconstruction,

IEEE

IMAGE

PI

Frieden,

B. R., Estimation

RECONSTRUCTION

- A New Role for Maximum (1977).

Entropy,

in SPSE

Con-

ference Proc. R. Shaw, ed., pp. 261-265

PI

Kikuchi, entropy

PI

Minerbo, G., MENT : A maximum entropy algorithm for reconstructing a source from projection data, Computer Graphics and Image Processing10, 48-68 (1979).

R. and Soffer, B. H., Maximum entropy expression, J. Opt. Sot. Am. 67, 16561665

Sanderson, PO1methods,

PI

939952

I. The

of fuel pin bundles by a maximum , 2685-2686 (1979).

entropy

IEEE Trans. Nuclear Science NS-26

Shepp, L. A. and Vardi, Y., Maximum likelihood reconstruction tomography, IEEE Trans. Medical Imaging MI-l, 113-122 (1982).

1121Jaynes, PI

J. G., Reconstruction

image restoration. (1977).

E. T., On the rationale (1982).

of maximum-entropy

Gull, S. F. and Skilling, J., Maximum Proc. IS,646659 (1984).

B. 1141Frieden, tomographic

R. and Zoltani, reconstruction,

X., [151 toZhuang? maximum

Ostevold, entropy

entropy

method

methods,

for emission

IEEE Proc. 70,

in image processing,

C. K., Maximum bounded Applied Optics 24, 201-207

entropy: (1985).

application

IEEE to

E. , and Haralick, R. M. , A differential equation approach image reconstruction, IEEE Trans. ASSP 35, 208-218

(1987).

I161 Park,

K. H. and Park. S. B., Maximum entropy image reconstruction for an object with opaque obstructions, IEEE Trans. Medical Imaging MI-6, 308-312 (1987).

1171Luenberger,

D. G., Optimization

by Vector

Space Method,

(Wiley,

New York,

1969).

PI

Dines, K. A. and Kak, A. C., Ultrasonic Ultrasonic Imaging, 1, 16-33 (1979)

83

attenuation

tomography

of soft tissue,