Multicriterion image reconstruction and implementation

COMPUTER VISION, GRAPHICS, AND IMAGE PROCESSING 46, 131--135 (1989)

Multicriterion Image Reconstruction and Implementation YUAN MEI WANG

Department of Scientific Instrument Engineering, Zhejiang University, People's Republic of China AND

WEI XUE LI) Zhejiang University, Hangzhou, People's Republic of China Received January 5, 1987; revised August 30, 1988 In this paper, a multiobjective optimization method of the m a x i m u m entropy image reconstruction from projection is described. We apply a new iterative algorithm to solve this problem. Computer simulation results are given. 9 1989 Academic Press, Inc. 1. I N T R O D U C T I O N

The problem of reconstructing an image from projections has arisen in a large number of scientific and engineering fields. G. T. Herman et al. [1-3] have majored in a comprehensive study of quadratic optimization methods for image reconstruction and obtained many results. The maximum entropy image reconstruction was formulated as the solution of constrained optimization problem [4-5]. The author [6-8] has suggested that maximum entropy image reconstruction from projections was represented as the solution of the nonlinear goal programming problem. In this paper, we discuss a new multiobjective model and iterated algorithm for image reconstruction in detail. Furthermore, the related property of goal programming was described. 2. N E W M O D E L A N D A L G O R I T H M

In the previous works [4-6], the maximum entropy image reconstruction from projections was formulated as the solution of the constrained optimization problem:

maximizeH(f) = f(x, y)>O

f f J ( x , y)ln(f(x, y))-i dxdy

over the set of continuous functions

f(x, y),

g(r, O) = f f J ( x , y)8(x

subject to the linear constraints:

cos 0 + y sin0

= f / ( x , y) ds

(1)

- r) dxdy (2)

131 0734-189X/89 $3.00 Copyright 9 1989 by Academic Press, Inc. All rights of reproduction in any form reserved.

132

WANG AND LU

and f ( x , y ) > O, for f ( x , y ) ~ D, where the line L = {(x, y ) : x c o s 0 + y s i n 0 = r} ( r , O ) ~ D a- {(r, 0 ) : - o o

< r <

oo, O < O < ~r)

8(.) is the Dirac delta function. 2.1. Multiobjective Optimization Model for Image Reconstruction In the multiobjective optimization model of the maximum entropy reconstruction, the resultant problem formulation is then given by the discrete version: Find f = (fl, f 2 , . . . , fmO T SO as to minimize a = [ ( P l ) , ( P : I + n21) . . . . . (P2,,, 2 + n2m2),(P3 + n3)]

(3)

such that fTln f + nl - Pl = - I n ( m 2 ) m 2

gi ---- E

q i j f j -F n 2i --

P2i = O,

j=l

i = 1,..., m 2

(4) (5)

IIg - Qfl[2 + n3 _ p3 = 0 and f , n = (nl, n21 . . . . , n2m2, n3) T,

p = (pl, p21,..., p2m2, p3) T >---0,

(6)

where g, f , Q are the known projection vector, image vector, and projection matrix, respectively. 2.2. The Solution Set The general form of image reconstruction represented as the multiobjective programming is: Find f = (fl, f2 . . . . . fN) T so as to minimize a=

{al(n,p)

..... a2(n,p)

}

(7)

such that hi(f)

+ ni - Pi = bi,

i = 1 ..... k

(8)

and f , n, p > 0,

(9)

where Eq. (7) is our achievement function and Eqs. (8) are the problem objectives. A decision variable, denoted as fj (with j = 1, 2 . . . . . N) will be assumed nonnegative unless otherwise noted. A deviation variable reflects the underachievement (negative deviation and denoted as ni, i = 1, 2 , . . . , k) or overachievement (positive deviation

MULTICRITERION IMAGE

133

a n d d e n o t e d as p~, i = 1, 2 . . . . . k ) of an objective i. A l l d e v i a t i o n variables are a s s u m e d to b e normegative unless otherwise specified. P r i o r to p r e s e n t i n g o u r discussion o f the n e w iterative a l g o r i t h m for solving the n o n l i n e a r m u l t i o b j e c t i v e p r o g r a m m i n g , we shall first describe s o m e of the definit i o n s a n d concepts. DEFINITION 1 (feasible solution). t u t e a feasible solution.

A n y set o f n o n n e g a t i v e fj, ni, Pt values consti-

DEFINITION 2 ( i m p l e m e n t a b l e solution). A n i m p l e m e n t a b l e s o l u t i o n is a feasible s o l u t i o n in w h i c h all a b s o l u t e objectives are satisfied. DEFImTION 3 (achievement f u n c t i o n o r vector). T h e goal p r o g r a m m i n g achievem e n t f u n c t i o n ( a ) is an o r d e r e d vector o f a d i m e n t i o n equal to the n u m b e r ( k ) o f p r e e m p t i v e priorities within the p r o b l e m a n d expresses the level o f achievement of e a c h o b j e c t i v e set within a given priority. DEFINITION 4 ( o p t i m a l solution). T h e s o l u t i o n ( f ) to a given goal p r o g r a m m i n g m o d e l is c o n s i d e r e d o p t i m a l if, for this s o l u t i o n ( t e r m e d f * ) , the c o r r e s p o n d i n g v a l u e o f a ( t e r m e d a*) is the s a m e or p r e f e r r e d to the value o f a for a n y other f e a s i b l e solution. 2.3. Algorithm W e h a v e e m p l o y e d a set of n e w iterative a l g o r i t h m for solving i m a g e reconstruct i o n p r o b l e m s w i t h highly satisfactory results. O u r iterative a l g o r i t h m is o n e of a class o f a c c e l e r a t e d search m e t h o d s . Such a l g o r i t h m s increase their search step size i f p r e v i o u s searches have b e e n successful a n d m a i n t a i n o r decrease the step size otherwise. W e m a y specify the steps of o u r iterative a l g o r i t h m for i m a g e r e c o n s t r u c t i o n as follows: NEW ITERATIVE A L G O R I T H M OF IMAGE R E C O N S T R U C T I O N Notation: f i m a g e vector; g p r o j e c t i o n vector; Q p r o j e c t i o n matrix; fr r e c o n s t r u c t e d image; N = m2 n u m b e r of projections. (1) Begin set k = 0. 1 = 0, y = h0, h ~ [0,1], 0: the primary set of increments; INm~x: the maximum iterative number; 0mln: the minimum of 0; 8: the minimum resolution of achievement function a. f(1) = f(0) = the arithmetic mean of all the rows of the matrix Q. Obtain the achievement function value a(f ~1)) for f~). Set ~a,0 =/(~)(2) Set j = k + l , k = k + l , l = l + l. (3) Examine the changes concerning ~j,0 so as to determine ~j, i as follows (a) set i = 1 Co) If a(~j, 0 + yi) < a(~j,o) then ~j,i = ~j,0 + Yi and go to 3(d) else go to 3(c). (c) If a(~j, 0 - yi) < a(l~j,o) then ~j,i = ~j,0 - Yi and go to 3(d) else ~j,i = ~j,0 and go to 3(d). (d) If 1 = m 2, go to 4 else i = i + 1 and go to 3CO). (4) If a(~j,m=) + 8 < a(f (j)) then f(p+l) = ~j.,,2 and go to 5 else f I N ~ then go to 8 else ~j.0 = (2x)f ~j) =f(j-1), h ~ [O,1] and go tO 2. (6) Set y = 2-xr, h ~ [0,1]. (7) If y < 0min then go to 8 else j = p + 1 and ~j,0 = f ( p + l ) and go to 2. (8) End of algorithm; select the last base point as the reconstructed results.

134

WANG AND LO TABLE 1 Algorithms

6~

Kashyap and Mittal [3] Least squares [3] Summation [3] Improved [3] Our GP model and algorithm

0.2445 (7) 0.2462 (7) 0.2738 (6) 0.2268 (8) 0.2012 (6)

3. COMPUTER SIMULATION RESULTS

The effectiveness of the new model and algorithm was shown by numerical examples of reconstruction for Shepp-Logan head data. A comparison to single objective optimization methods were given. The advantages of our modeling and algorithm for image reconstruction are manifesting themselves in combination with other optimization methods. The measures of image quality defined below in the reconstruction literature (3): N

~

Z

}1/2

/~j~l

N

Z (fo(i, j l - f,(i, J))"

i=1 j = l

(fo(i, j ) - fo) 2 i=

,

'=

where fo(i, j ) and f,(i, j) denote the gray levels in the ith pixel of the j t h row of the digital 256 • 256 test phaton and the reconstruction, respectively, and f0 denotes the average gray level in the digital test picture. The simulation results are reported in Table 1 and in Fig. 1, respectively. For each the algorithm results of the sixth iteration is reported. 4. CONCLUSIONS

The multiobjective programming method for image reconstruction has many advantages over single quadratic optimization in terms of storage, algorithm simplicity, and convergence rate. We add an entropy constraint to the classical image reconstruction in quadratic optimization model in order to preserve a desired level of accessibility in the solutions to the image reconstruction. We add an entropy of the image of the objective function with the purpose of introducing an element of

(a)

(b) FIG. 1. Phantom and multicriterion reconstruction.

MULTICRITERION IMAGE

135

smoothing into the solution of nonlinear programs; the element of the optimal solutions become strictly positive. The entropy constraint may in such cases be looked upon as a substitute for lost complexity. Our modeling and algorithm may be applied to the sonar and SAR multidimensional signal reconstruction, and multidimensional power spectrum estimations. REFERENCES 1. G. T. Herman and A. Lent, Quadratic optimization for image reconstruction, I, Comput. Graphics Image Process. 5, (1976), 319-332. 2. E. Artzy and G. T. Herman, Investigation of quadratic optimization techniques for image reconstruction, in Proceedings, 1977 IEEE Conferehce on Decision and Control, New Orleans, Louisiana, Dec. 1977, pp. 350-360. 3. E. Artzy, T. Eliving, and G. T. Herman, Quadratic optimization for image reconstruction, II, Comput. Graphics Image Process. 11, 1979, 242-261. 4. S. J. Wernecke et al., Maximum entropy image reconstruction, IEEE Trans. Comput. C-26, 1977, 351-364. 5. G. Minerbo, MENT: A maximum entropy algorithm for reconstructing a source from projection data, Comput. Graphics Image Process. 10, 1979, 48-68. 6. Y. M. Wang, Goal Programming method of maximum entropy image reconstruction from projections, presented at International Conference on Numerical Optimization and Applications, Xian, Shaanxi, China, June 1986. 7. Y. M. Wang, Nonlinear goal programming model and algorithm of maximum entropy image reconstruction from projections, presented at the International Coeference on Acoustics, Speech, and Signal Processing, Beijing, Institute of Acoustics, Academia Sinica, 1986. 8. Y. M. Wang, New method of maximum entropy image reconstruction, presented at the International Conference on Acoustic, Speech, and Signal Processing, Tokyo, Japan, April 1986.