Two dimensional phase retrieval using neural networks

T W O DIMENSIONAL PHASE RETRIEVAL USING NEURAL NETWORKS Adrian Burian, Pauli Kuosmanen, Jukka Saarinen, Corneliu Rusu Digital Media Institute, Tampere University of Technology PO Box 553, FIN-33101, Tampere, Finland E-mail: burian, pqo, jukkas,cornelin @cs.tut.fi

Abstract. The two-dimensional (2-D) phase retrieval problem is to reconstruct an image from only the spectrum magnitude. This problem emerges when the phase of the 2-D signal is apparently lost or impractical to measure. For 2-D spatially limited non-negative objects, characterized by analytic spectrum, the solution is unique. In this paper we propose the use of neural networks for solving the 2-D phase retrieval problem. The neural network incorporates a combination of the maximum entropy estimation algorithm with

some additional nonlinear constraints. These constraints make use of additional unknowns related to the real and imaginary parts of the image spectrum. The solution is provided by the steady state of the neural network, then is verified and eventually improved with the iterative Fourier transform algorithm. The,obtained simulation results demonstrate the high efficiency of the proposed approach.

INTRODUCTION Phase information has a fundamental importance in many two-dimensional (2-D) signal processing problems. Different applications, such as electron microscopy, magnetic resonance imaging, wavefront sensing, optical imaging, requires the reconstruction of an image known to have a compact support, from its Fourier transform magnitude. The image is reconstructed if the missing Fourier phase is recovered - hence the name phase retrieval [6]. In the case of 1-D signals the phase retrieval problem has an unique solution only for minimum-phase signals, otherwise the problem is 'ill-posed' [2]. For 2-D signals it is known that the solution is unique, but does not exist yet a precise algorithm to derive it. Because of its importance, the phase retrieval problem has attracted many research efforts, and a number of algorithms have been proposed in the literature. The most common one is to use an iterative Fourier transform (IFT) algorithm, which alternates between spatial and transform domains [5]. However, these algorithms usually stagnate, failing to converge to a solution. In [15] a 'battery' of digital signal processing 0-7803-6278-0/00/$10.00 (C)2000 IEEE

652

algorithms was used to produce an algorithm that comes close to expressing the solution in closed form, while in [13] a window function was used. In [l]a maximum entropy method was considered, and the solution of the considered phase problem was obtained numerically, by using steepest-descent method. Usually, the used methods are combined with the IFT algorithm in order to improve the probability of success of phase retrieval [9]. In this paper we consider 2-D phase retrieval using neural networks. The maximum entropy algorithm with additional nonlinear constraints is used in this purpose. Maximum entropy methods are being applied to a range of problems involving incomplete information. Examples of such applications include spectral estimation algorithms and algorithms for image reconstruction. The essence of the maximum entropy method is to maximize the entropy of the reconstructed distribution, subject to satisfying constraints on the distribution. In the literature there exist several proposals for architectures of neural networks for maximum entropy estimation. In [4] the entropy was maximized by a constrained gradient ascent algorithm that was implemented by an artificial neural network architecture. This architecture contained three layers: the probability layer, the gradient layer, and the constrained layer. Each unit from the gradient layer receives input from a single unit from the probability layer. The constraint units receive input from the gradient units, and pass their outputs to the probability units. The probability units compute the sum of the products of the states of the constraint units and the weights on the corresponding connections. This sum is then added to the current state. The network was allowed to run cyclically until the states of the units stabilized. The authors used this architecture for prognosis of breast cancer patients. In [7] a neural network implementation of the maximum entropy reconstruction algorithm was proposed. The degraded image was assumed to be created by a process of convolution of the original image with a blurring operator, and the addition of Gaussian noise. With the advances in new technologies, and especially in very large scale integration - VLSI - technology, the use of neural networks in solving optimization problems has been proposed. There exist also proposed circuit structures that can be used, after slight modifications, in related problems [3]. The neural network approach enables us to solve these problems in real time, due to the massively parallel operations of the computing units and faster convergence properties. We have incorporated the maximum entropy estimation and the respective nonlinear constraints into the neural network structure. The advantages of using a derivation of the maximum entropy algorithm consist in good extrapolation features and high stability to noise. The disadvantage of the maximum entropy method is the massive computational requirements. We overcome this disadvantage by using a highly parallel neural network computational structure.

0-7803-6278-0/00/$10.00( C )2000 IEEE

653

MAXIMUM ENTROPY AND 2-DPHASE RETRIEVAL Let the required reconstructed image has pixel values represented by positive numbers x1 ,. . . ,x n which are to be determined, and on which the entropy is given by:

where p i denotes the probability of encounter of a single photon with pixel i. For the phase retrieval problem, the energy of the image is constant during the reconstruction process, and a linear relation between S(p) and S(x) exist. The intensity of a pixel is proportional to this probability. This is one way to define the entropy. Another one uses both bounds of the intensities of the original image. If these are known: 0 5 xi 5 T , the following expression for the entropy can be derived: n

s ( P ~ ,,*P*n .>= - C b i h ~+i ( T - P ~ I M T - P ~ I I -

(2)

i

In the sequel, the 2-D phase retrieval problem will be discussed. The 1 ) of , size Fourier modulus IX0(nlI C ) / of a real and nonnegative object ~ ( m N x N is considered to be the given data. The reconstruction problem is to find a real and nonnegative solution such that I X ( n ,k)l = IXo(n,k ) l . The 2-D unitary discrete Fourier transform pair is defined as:

1.

n=O

k=O

Given the modulus IX(n, IC)! = IXo(n, k)I, the reconstruction of the image is equivalent to the reconstruction of the Fourier phase 4(nl IC) from the magnitude. The phase retrieval problem is equivalent with the autocorrelation retrieval: &(n, k) = IDFT [IXo(~,~)l21.

(5)

In order to solve the problem, the usual approach is to find an iterative solution z(m,Z) such that I X ( n ,k)l- > IXo(n, I C ) ! . An important question that is needed to be raised is that whether the obtained solution with I X ( n , k ) l - > IXo(n,k)l will lead to x(m,I)- > xo(m,Z). Unfortunately, a direct and complete answer to this question does not exist yet [9].Instead, the question of the uniqueness of the solution of 2-D phase retrieval problem 0-7803-6278-0/00/$10.00 (C)2000 IEEE

654

was thoroughly investigated. It has been established that the solution of a 2-D phase retrieval problem, if it exist, is typically unique to within a few trivial ambiguities - change of sign, shift in position, reversal of orientation of the solution. Previous reported results shown that the 2-Dphase retrieval problems are usually solvable and unique [ll],[13], [15]. We separate the real and imaginary parts on both sides of equation (3). Let us denote by Mnk the known spectral magnitudes, and by a$ and b$ the constants related to discrete Fourier transform: a nm[ k

- cos(2n(mn + l k ) / N ) / N ,

b$ = - sin(2n(mn

+Ik)/N)/N.

(6) (7)

The unknown phase is denoted by &k. The unknown nonnegative object x ( m ,1 ) is denoted by x m [ . The following equations are obtained:

Mnk Sin &k

xmrbkt.

= m

(9)

l

The magnitudes of the spectrum being known, also the sum of all values to be reconstructed is known CmClxm1 = Moo. Let us consider the following nonnegative variables: tnk = 1 Snk

+

=1

(10)

COS$&,

+ sin &k.

(11)

By using the above definitions, the phase retrieval problem can be written as a nonlinear optimization problem, by the way of the maximum entropy restoration: find max

-x m

ln(tnk) + Snk h(snk)],

X m l h(Xm1) -

l

n

(12)

k

subject to a) Zml,tnk,Snk 3 0;

0-7803-6278-0/00/$10.00 (C)2000 IEEE

655

(13)

with normalization:

By analyzing the above relations, we can observe that the complementary problem can also be defined in the same way. In this case the phase represent the given data and the spectral magnitudes are unknown. The phase normalization condition (16) is not needed anymore and the maximum entropy expression (12) is modified accordingly (Mnk instead of Snk and tnk). In order to improve the probability of success of the phase retrieval, the IFT algorithm is applied, with the obtained image as starting point. To use only the IFT algorithm is not a good idea, because of the unpredictability of its convergence. In [lo] a potential source of error in the numerical implementation of the algorithm is discussed. In [ll] it was noted that the IFT algorithm might stagnate at some ambiguous image, which may be similar or distinctly dissimilar in appearance to the true object.

NEURAL NETWORK BASED 2-DPHASE RETRIEVAL The maximum entropy method is well described in the literature, but the standard algorithm has the disadvantage of computational inefficiency. Also, the 2-D phase retrieval existing algorithms are very complex. Since an artificial neural network has a strong computational capability, it can be used to solve these problems. As the recurrent neural networks and the entropy are highly nonlinear, the energy minimization procedure should be modified in order to get the maximum entropy approach [12], [14]. The nonlinear optimization problem (12) can be written in a matrix form as follows: minimize f(y),

(17)

subject to the equality constraints:

QY = C ;

(18)

and the simple bound constraints:

with normalization: (Ynk

- 1)' t (ynk+pp

- 1)2 = 1 for

N 2 < n, k 5 2 N 2 .

In the above relations, we used the following notations:

. f ( ~= ) -s(Y) =

Cz/jh y j , Yj 2 0,vji j

0-7803-6278-0/00/$10.00 (C) 2000 IEEE

656

(20)

and

Q E R2N2x3Nz, Y E R 3 N ZC, E R2N2. The matrices A E R N Z x Nand Z B E R N Z X Ncontains 2 the constants given by relations (6) and (7) respectively, arranged on lines after nk and on columns after ml. The column vector M E RNZcontains all the known magnitudes. For example, for 256 gray level images, the superior bound for the corresponding X values is y3 max = 255 for j = 0 , 1 , . ,NxN - 1 and for the corresponding T and S values is y j = 2 for j = NxN, . ,3xNxN - 1. These box constraints are realized by employing limiting integrators with nonlinear limiters at their outputs. In this paper we will solve the 2-D phase retrieval problem by incorporating the above relations into the structure of a neural network. If we consider a neural network with fully interconnected neurons, an energy function U and an entropy function S can be associated with it. For a positive convex entropy function, and for any non-decreasing temperature sequence T, the neural network admits a Lyapunov function, which is the Helmholtz free energy of the system [7], [8]:

F=U-TS.

(24)

The neural network will evolve in time until it reaches an equilibrium state that corresponds to a minimum of the free energy function F, which simultaneously minimizes the energy and maximizes the entropy. The additional constraints were incorporated into the neural network by the way of the energy function, which was modified in order to suit the available information. The nonlinear, sigmoidal transfer function that determines the relation between an input vi and an output yi is given by: yi = 1

+ tanh

(7)

The slope of the transfer function at the inflection point vi = 00 (where vo is the offset) constitute the maximum gain of the amplifier in the practical realization, and is given by:

0-7803-6278-0/00/$10.00 (C)2000 IEEE

657

Figure 1: Example of image reconstruction of a 8 x 8 pixel resolution text letter: a) original image, b) the reconstructed image after 50 iterations, c) the obtained solution after 100 iterations, d) the shifted reconstructed image. In order to obtain the phase normalization, the relation (25) is modified as follows:

%T,S

=1

+ tad2

(

( V~T;VQ)

viT9:-vo

)

,

+ tanh2 ( ~ p v ~ )

(27)

for N2 5 iT < 2N2, 2N2 5 is < 3N2. By inserting the phase normalization condition in this manner, the neural network still admits a Lyapunov function. This is because the function involved is positive monotone increasing and differentiable. We encoded the constraints into function U in the following way:

where kl is a scaling parameter. We obtained:

The time evolution of the neural network is described by:

The steady state of the neural network with the nonlinear constraints satisfied, will provide the solution for the 2-D phase retrieval problem. We considered that the constraints are hard, i.e., we required that they are satisfied as closely as possible, in order to guarantee a minimum error, and to avoid the problems with the maximum entropy solutions with soft constraints. The obtained solution is verified and eventually improved with the IFT algorithm. The simulations we made demonstrate that this procedure ensures reliable convergence to the solution.

0-7803-6278-0/00/$10.00 (C) 2000 BEE

658

Figure 2: Example of image reconstruction for a 32 x 32 pixel resolution image: a) original image, b) the reconstructed image after 500 iterations, c) the obtained solution after 1000 iterations.

SIMULATION RESULTS To illustrate the above approach, the first used image was a synthetic one, of a 8 x 8 pixel resolution text letter shown in Figure la). The magnitudes of the 2-D Fourier transform of this image are the given data. Figure l b ) shows the reconstructed image after 50 iterations, and Figure 2c) the final reconstructed image. It is obvious that this image is the shifted version of the original one. Only one temperature was used during the whole process T = 0.01. We repeated the same procedure using a 32 x 32 part of the Lena image. Simulation results for this image are shown in Figure 2. In Figure 2.a) the original image was presented. The reconstructed image after 500 iterations is presented in Figure l.b), and after 1000 iterations in Figure 1.c). Starting temperature was T = lo-*. For the presented simulation it was raised in 4 steps until T = 1. The raising of the temperature was decided by a very small decreasing value of the error during 50 iterations. The small size of the used images is justified by the high dimensionality of the involved matrices. But the obtained results are also valid for bigger sizes, which are not a problem for the hardware implementations. The power of the neural network implementation lies in its high parallelism. For the used neural network, it wi~simportant that the data presented to his nodes to be scaled appropriately. The entire neural network’s output data were scaled between 0 and 2. This was done in order to encourage equitable distribution of importance. The variables with large absolute magnitudes tend to exercise more influence on the network’s response than those having small magnitudes. When the outputs corresponding to the magnitudes of the reconstructed image were considered integer numbers between 0 and 255, the network had a decreasing error function only when different temperature parameters, for different output scales, were used. But this correction takes time, because there are more parameters to adjust and sometimes the convergence may fail to occur, Also, at each iteration, the data corresponding to the magnitudes were scaled on 256 levels between 0 and 2. Only one final scaling on these levels it is not enough, because it can cause a high increase of the error. 0-7803-6278-0/00/$10.00 (C) 2000 IEEE

659

The error of the phase retrieval problem was computed as the mean square differences between the obtained magnitude spectrum, and the given magnitude spectrum. In the simulations we made, the neural network stagnate sometimes at some bigger values of the error, but in all of these cases, after a small number of IFT algorithm iterations (about 20), the correct solution was obtained.

CONCLUSIONS In this paper, we have addressed the image reconstruction from its Fourier transform intensities by using neural networks. Reconstruction of an image from its Fourier magnitude is regarded as a nonlinear optimization problem. A derivation of the maximum entropy estimation algorithm with some nonlinear constraints for the real and imaginary parts of the spectrum data was included into the structure of a neural network. The intrinsic parallelism inherent in a neural network makes it possible to overcome the low speed and high number of computations of the maximum entropy method used for the 2-D phase retrieval problem.

REFERENCES [l] A. Bajkova, “Phase Retrieval Algorithm based on Maximum Entropy Method,” Proc. of Int. Symp. on Signals, Systems and Electronics, pp. 307-309, 1995. [2] A. Burian, P. Kuosmanen and C. RUSU,“1-D Non-Minimum Phase Retrieval by Gain Differences,” Proc. of IEEE Int. Conf. on Electronics, Circuits and Systems, pp. 1001-1005, Sept. 1999. [3] A. Cichocki and R. Unbehauen, Neural Networks for Optimization and Signal Processing, John Wiley & Sons, 1993. [4] C. desilva and P. Choong, “A Network Architecture for Maximum Entropy Estimation,” Proc. of IEEE Int. Conf. on Neural Networks, vol. 3, pp. 1438a-1442, 1994. [5] J. Fienup, “Phase retrieval algorithms: a comparison,” Applied Optics, vol. 21, no. 15, pp. 2758-2769, 1982. [6] M. Hayes, “The reconstruction of a multidimensional sequence from the phase or magnitude of its Fourier transform,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 2, pp. 140-154, 1982. [7] D. Ingman and Y. Merlis, “Maximum Entropy Signal Reconstruction with Neural Networks,” IEEE Trans. on Neural Networks, vol. 3, no. 2, pp. 195-201, March 1992. [8] G. Lendaris, K. Mathia and R. Saeks, “Linear Hopfield Networks and Constrained Optimization,” IEEE Trans. Systems, Man, and Cybernetics - Part B: Cybernetics, vol. 29, no. 1, pp. 114-118, Feb. 1999. (91 Z. Mou-yan and R. Unbehauen, “Methodsfor Reconstruction of 2-D Sequences from Fourier Transform Magnitude,” IEEE Trans. Image Processing, vol. 6, no. 2, pp. 222-233, Feb. 1997. 0-7803-6278-0/00/$10.00 ( C ) 2000 IEEE

660

[lo] J. Sauz, T. Huang and T. Wu, “A note on iterative Fourier transform

phase reconstruction from magnitude,” IEEE Trans. on Acoustic, Speech, and Signal Processing, vol. 32, no. 6, pp. 1251-1254, 1984. [ll] J. Seldin and J. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” Journal Opt. Soc. Amer., vol. 7, no. 1, pp. 3-13, 1990. [12] J. van der Berg, “The Most General namework of Continuous Hopfield Neural Network,’’ Proc. of Int. Work. on Neural Networks for Identification, Control, Robotics, and Signal/Image Processing, pp. 92-100, 1996. [13] K. Wooshik, ‘:Two-Dimensional Phase Retrieval using a Window Function,” Proc. of IEEd Int. Conf. on Acoust., Speech and Signal Processing, vol. 3, pp. 1625-1628, 1999. [14] Y. Xia, J. Wang and D.Hung, “Recurrent Neural Networks for Solving Linear Inequalities and Equations,’’ IEEE Trans. Circuits and Systems-I: Fundamental Theory and Applications, vol. 46, no. 4, pp. 452-462, April 1999. [15] A. Yagle, “Digital Signal Processing Solutions to 2-D Phase Retrieval Problems,” Proc. of IEEE Int. Conf. on Image Processing, pp. 712-716, 1998.

0-7803-6278-0/00/$10.00(C) 2000 IEEE

661