- Email: [email protected]
Relevance for blind deconvolution of recovering Fourier magnitude from phase
1
OPTICS
RELEVANCE FOR BLIND FROM PHASE
COMMUNICATIONS’
DECONVOLUTION
1 July 1987
OF RECOVERING
FOURIER
MAGNITUDE
R.G. LANE and R.H.T. BATES Electrical and Electronic Engineering Department, University of Canterbury, Christchurch I, New Zealand Received
28 October
1986; revised manuscript
received
10 March
1987
Blind deconvolution of a single blurred positive image is shown to be achievable with the aid of an algorithm for recovering Fourier magnitude from phase (modulo IT). The algorithm, which is illustrated with computational examples, is shown by example to be robust in the presence of noise. While the magnitude-from-phase (modulo 2x) problem is confirmed to be efficiently solvable using Fienup-type iterative algorithms, it is pointed out that invoking the concept of the zero-sheet of the Fourier spectrum complements Hayes’s previous argument that more constraints are needed to ensure uniqueness (e.g. the image support must be known a priori) than for the phase-from-magnitude problem.
It is widely held that recovery of the magnitude of a visibility (by which is here meant the Fourier transform of an image) from its phase is much less difficult than interferring the phase from the magnitude[ 1,2]. While this is true in one sense, it is false in another as we demonstrate below. We have been led to re-examining the recovery of magnitude from phase by studying the blind deconvolution of a single blurred image. The problem is: recover an unknown original imagefix) from a given blurred image g(x) without prior knowledge of the point spread function ( psf) h(x) , where g(x) =&)-h(x)
9
.
(1)
with x being the position vector of an arbitrary point in two-dimensional image space and ‘-’ being the convolution operator. It is here understood that both f(x) and h(x) are positive (i.e.real and non-negative) and have compact support (but neither support is given a priori). If u is the position vector of an arbitrary point in Fourier space, and upper-case letters identify the visibilities of images denoted by lower-case letters, then the convolution theorem [ 31 indiates that G(u)=F(u) We have employing
H(U).
(2)
found that, starting with 1G(u) 1 and Fienup’s phase recovery algorithms [ 41,
0 030-4018/87/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
using a series of initial pseudo-random phase distributions, it is always (i.e. in all the cases we have tried) possible to generate an estimate of G(U)=
eitherF*(u)H(u)
orF(u)H*(u)
,
(3)
where the asterisk denotes complex conjugation. It is important to recognise that we cannot be sure which version of G( U) has been reconstructed in any particular instance, because there is no a priori information available concerning the forms off(x) and h(x). It follows from eqs. (2) and (3) that dividing G(U) and c(u) gives twice the phase of either F(u) or H(U). Once a successful magnitude recovery scheme has been developed, both f(x) and h(x) can be individually reconstructed. We cannot be sure (in principle), which is which, but in practice it will usually be obvious. So it makes sense from now on to consider that we know 2 (phase{F( u) }), but with a vital proviso: since both F(u) and G(U) are generated computationally (as opposed to analytically), this doubled phase is only obtained modulo 2a. When it is divided by two, we have the quantity Y(u) =phase{F(u)}
modulo
n .
(4)
It is worth reiterating here the point made previously by Hayes [ 21 that, if either the image or the psf is conjugate symmetric, the equivalent phase inforB.V.
11
Volume 63, number
1
OPTICS
COMMUNICATIONS
mation is given by G(u)lG*(u) because the Fourier transform of a conjugate symmetric object is necessarily real. We conclude this paper by illustrating the recovery of 1F(u) 1 from Y(u), and we do so for noisy data. However, we first consider the general problem of recovering IF(u) I from phase {F(u)}. If Ax) has compact support, so that F(u) is an entire function, then F(u) is zero on a continuous surface, called the zero-sheet Z,, in the four-dimnensional space formed by generalising each of the two real components of u to a complex variable. While the necessity for this is argued in detail elsewhere [ $61, it is worth remarking here that a zerosheet can be formed by, first, computing the point zeros lying in any two-dimensional plane in the fourdimensional space and, second, following the tracks of these zeros through all such planes [ 61. The zerosheet must be continuous because any discontinuity is incompatible with F(u) being an entire function. Given only its zero-sheet, F(u) can be reconstructed to within a complex constant. It should be noted that the compactness of the support off(x) precludes the possibility of F(u) being multiplied by an entire function which goes to zero only at infinity in the aforementioned four-dimensional space (e.g. a gausSian). It is possible to compute Z, from )F(u) I when the latter is sampled sufficiently finely to compute the autocorrelation of_/(x) [ 561. When only phase {F(u)} is given, it is impossible to infer ZF uniquely, since it is always possible to multiply F(u) by an entirely positive visibility Q(u), e.g. p(u) =F(u)Q(u)
>
(5)
in which case phase{F(u)} =phase{F(u)}. Furthermore, the support ofnr) is larger than the support offlx) . The support of q(x) is effectively finite, since the given samples of Q( u) can be taken to satisfy the sampling theorem [ 31. It is then inevitable that Q(u) is characterised by a zero-sheet Z, existing in the finite part of the four-dimensional space [ $61. So, zpz,uz,
)
(6)
from which we deduce the perhaps surprising result that while f(x), or more precisely its image-form [ 7,8], is in general uniquely defined by IF(u) 1, we can only deduce the imageflx) uniquely from phase 12
1 July 1987
{F(u)} if we also know that the support off(n) is the smallest consistent with the given phase distribution, in which case q(x) reduces to a two-dimensional delta function. The magnitude-from-phase problem is thus “less unique” than the phase-from-magnitude problem. However, when the minimum support constraint is known to apply, algorithms of the Fienup-type[4] recover I F(u) I from phase {F(u)} much more effectively than vice versa, and they can handle complexvalued images as efficiently as positive images (which is not true for the phase-from-magnitude problem). The phase must be over-sampled by comparison with what is necessary when both I F(u) I and phase {F(u)} are given, just as the magnitude must be oversampled for the phase-from-magnitude problem, in order that the reconstructed image can be forced to zero outside the minimum support. We feel that the concept of the zero-sheet provides an illuminating geometrical interpretation of the magniproblem, tude(phase)-from-phase( magnitude) which nicely complements Hayes’s detailed algebraic arguments [ 21. From a pragmatic algorithmic viewpoint, the main contribution of this paper is the strategy (described in detail below) for recovering the Fourier magnitude from the modulo II Fourier phase. When the image support is unknown a priori, but the minimum support constraint is known to apply, the Fienup-type algorithms are applied with successively decreasing image-boxes (an image-box is that rectangular region of image-space, centred at the origin with sides parallel to arbitrarily chosen Cartesian axes, outside of which the image is set to zero) until a large increase in the pixel level outside the support is observed. This indicates that the support used previously is the smallest support consistent with the data. Given w(u), as defined by (4), rather than phase{F(u)} itself, we have found the following approach to reconstucting IF(u) I to be effective. The iterations are started by computing a pseudo-random positive image e(x) which vanishes outside a chosen image-box b(x). We next compute the visibility of e(r), which is conveniently written as A(u)exp[i@u)J, where A(u) is positive and O(u) is real. O(u) is then subtracted (modulo2lr) from w(u). At each pixel for which I v(u) -13(u) I is less or
Volume 63, number
1
Fig. 1. Simple image consisting
OPTICS
COMMUNICATIONS
of five isolated pixels.
greater, respectively, than n/2 we redefine 19(u) as w(u) or (w(u) + n) . After recomputing an updated e(x) by Fourier transformation, we set to zero its negative pixels inside, and all of its pixels outside, b(x). An iterative loop is set up thereby. We assess the rate of convergence by calculating, at each iteration, the fraction e of ]e(r) I for which e(x) lies outside b(x) and is negative inside b(r). This procedure is repeated for image-boxes of decreasing size, until the condition noted in the previous paragraph is observed. Fig. 1 shows a simple image consisting of five isolated pixels. Figs. 2 and 3 show the images reconstructed when the linear dimensions of b(x) are set 2OWlarger and smaller, respectively, than those of the correct image-box, which is the b(x) just enclosing the f(r) shown in fig. 1. Fig. 4 shows the variation of 6 with N (the number of iterations) for these image-boxes and the correct image-box. Fig. 2 shows the blurring of f(x) which is inevitable when the image-box is larger than the correct b(x) - note that each pixel is similarly blurred, i.e. it is convolved with a q(x) as defined above. The curve of E versus N for the largest b(x) confirms that convergence is achieved, just as satisfactorily as for the correct b(x). We note that the reconstruction is effectively identical to fig. 1 when b(x) equals the correct image-box. When b(x) is smaller than the correct b(x), fig. 4 indicates that convergence is never achieved, as is also suggested by the sporadic appearance of fig. 3.
1 July 1987
Fig, 2. Reconstruction of image shown in fig. 1 when the linear dimensions of b(x) are 20% larger than those of correct b(x).
These results, together with other computational experience persuade us that we can estimate the linear dimensions of the correct b(x) from inspection of the curves of t versus N. Fig. 5 shows an image containing intricate detail. We added to the oversampled Fourier transform uncorrelated pseudo-random noise, uniformly distributed in the range [ 0, 2n] in phase and [ 0, m] in magnitude where m was 5Okof the DC level. The noise E(u) was forced to be conjugate-symmetric, implying that E( - U) =E* ( U) , thereby ensuring that
Fig. 3. Reconstruction of image shown in fig. 1 when the linear dimensions of b(x) are 20% smaller than those of correct b(x).
13
Volume 63, number
1
OPTICS
the image noise was real. The magnitude of the Fourier transform was discarded and its modulo II phase was retained. Fig. 6 shows the image reconstructed after 500 iterations (when t has decreased to 1.1 x 1Oe3) when b(x) was made equal to the correct image-box. Given the significant noise level and the consequent degradation of high spatial frequency information, fig. 6 suggests that Fienup-type algorithms for recovering magnitude from phase (modulo n) are practical and robust in the presence of noise. This result is comparable to those obtained
Fig. 6. Reconstruction (sampled twice as finely as demanded by the sampling theorem in both directions in Fourier space) to which noise had been added as prescribed in final paragraph of text.
with the conventional Fienup algorithms when applied to the recovery of Fourier phase distributions from noisy magnitude data [ 8,9].
References
Fig. 5.64 x 64 pixel positive image, quantised levels from black to white.
14
to 256 equi-spaced
[ I] A.V. Oppenheim and J.S. Lim, Proc. IEEE 69 (198 1) 529. [2] M.H.Hayes, IEEETrans. ASSP-30 (1982)140. [3] R.H.T. Bates and M.J. McDonnell, Image restoration and reconstruction (Clarendon Press, Oxford, 1986). [4] J.R. Fienup, Appl. Optics 21 (1982) 2758. [S] R.G. Lane, W.R. Fright and R.H.T. Bates, Direct phase retrieval, to appear in IEEE Trans. ASSP-35 (April 1987). [ 61 R.G. Lane and R.H.T. Bates, Automatic multi-dimensional deconvolution, to appear in J.Opt. Sot. Am. A. 4 (January (987). [7] R.H.T. Bates and W.R. Fright, J. Opt. Sot. Am. 73 (1983) 358. [ 81 R.H.T. Bates and D. Mnyama, in: Advances in electronics and electron physics, Vol. 67, ed. P.W. Hawkes (Academic Press, 1986) p. 1. [9] J.R. Fienup, in: Indirect imaging, ed. J.A. Roberts (Cambridge University Press, 1984) p. 99.
OPTICS
RELEVANCE FOR BLIND FROM PHASE
COMMUNICATIONS’
DECONVOLUTION
1 July 1987
OF RECOVERING
FOURIER
MAGNITUDE
R.G. LANE and R.H.T. BATES Electrical and Electronic Engineering Department, University of Canterbury, Christchurch I, New Zealand Received
28 October
1986; revised manuscript
received
10 March
1987
Blind deconvolution of a single blurred positive image is shown to be achievable with the aid of an algorithm for recovering Fourier magnitude from phase (modulo IT). The algorithm, which is illustrated with computational examples, is shown by example to be robust in the presence of noise. While the magnitude-from-phase (modulo 2x) problem is confirmed to be efficiently solvable using Fienup-type iterative algorithms, it is pointed out that invoking the concept of the zero-sheet of the Fourier spectrum complements Hayes’s previous argument that more constraints are needed to ensure uniqueness (e.g. the image support must be known a priori) than for the phase-from-magnitude problem.
It is widely held that recovery of the magnitude of a visibility (by which is here meant the Fourier transform of an image) from its phase is much less difficult than interferring the phase from the magnitude[ 1,2]. While this is true in one sense, it is false in another as we demonstrate below. We have been led to re-examining the recovery of magnitude from phase by studying the blind deconvolution of a single blurred image. The problem is: recover an unknown original imagefix) from a given blurred image g(x) without prior knowledge of the point spread function ( psf) h(x) , where g(x) =&)-h(x)
9
.
(1)
with x being the position vector of an arbitrary point in two-dimensional image space and ‘-’ being the convolution operator. It is here understood that both f(x) and h(x) are positive (i.e.real and non-negative) and have compact support (but neither support is given a priori). If u is the position vector of an arbitrary point in Fourier space, and upper-case letters identify the visibilities of images denoted by lower-case letters, then the convolution theorem [ 31 indiates that G(u)=F(u) We have employing
H(U).
(2)
found that, starting with 1G(u) 1 and Fienup’s phase recovery algorithms [ 41,
0 030-4018/87/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
using a series of initial pseudo-random phase distributions, it is always (i.e. in all the cases we have tried) possible to generate an estimate of G(U)=
eitherF*(u)H(u)
orF(u)H*(u)
,
(3)
where the asterisk denotes complex conjugation. It is important to recognise that we cannot be sure which version of G( U) has been reconstructed in any particular instance, because there is no a priori information available concerning the forms off(x) and h(x). It follows from eqs. (2) and (3) that dividing G(U) and c(u) gives twice the phase of either F(u) or H(U). Once a successful magnitude recovery scheme has been developed, both f(x) and h(x) can be individually reconstructed. We cannot be sure (in principle), which is which, but in practice it will usually be obvious. So it makes sense from now on to consider that we know 2 (phase{F( u) }), but with a vital proviso: since both F(u) and G(U) are generated computationally (as opposed to analytically), this doubled phase is only obtained modulo 2a. When it is divided by two, we have the quantity Y(u) =phase{F(u)}
modulo
n .
(4)
It is worth reiterating here the point made previously by Hayes [ 21 that, if either the image or the psf is conjugate symmetric, the equivalent phase inforB.V.
11
Volume 63, number
1
OPTICS
COMMUNICATIONS
mation is given by G(u)lG*(u) because the Fourier transform of a conjugate symmetric object is necessarily real. We conclude this paper by illustrating the recovery of 1F(u) 1 from Y(u), and we do so for noisy data. However, we first consider the general problem of recovering IF(u) I from phase {F(u)}. If Ax) has compact support, so that F(u) is an entire function, then F(u) is zero on a continuous surface, called the zero-sheet Z,, in the four-dimnensional space formed by generalising each of the two real components of u to a complex variable. While the necessity for this is argued in detail elsewhere [ $61, it is worth remarking here that a zerosheet can be formed by, first, computing the point zeros lying in any two-dimensional plane in the fourdimensional space and, second, following the tracks of these zeros through all such planes [ 61. The zerosheet must be continuous because any discontinuity is incompatible with F(u) being an entire function. Given only its zero-sheet, F(u) can be reconstructed to within a complex constant. It should be noted that the compactness of the support off(x) precludes the possibility of F(u) being multiplied by an entire function which goes to zero only at infinity in the aforementioned four-dimensional space (e.g. a gausSian). It is possible to compute Z, from )F(u) I when the latter is sampled sufficiently finely to compute the autocorrelation of_/(x) [ 561. When only phase {F(u)} is given, it is impossible to infer ZF uniquely, since it is always possible to multiply F(u) by an entirely positive visibility Q(u), e.g. p(u) =F(u)Q(u)
>
(5)
in which case phase{F(u)} =phase{F(u)}. Furthermore, the support ofnr) is larger than the support offlx) . The support of q(x) is effectively finite, since the given samples of Q( u) can be taken to satisfy the sampling theorem [ 31. It is then inevitable that Q(u) is characterised by a zero-sheet Z, existing in the finite part of the four-dimensional space [ $61. So, zpz,uz,
)
(6)
from which we deduce the perhaps surprising result that while f(x), or more precisely its image-form [ 7,8], is in general uniquely defined by IF(u) 1, we can only deduce the imageflx) uniquely from phase 12
1 July 1987
{F(u)} if we also know that the support off(n) is the smallest consistent with the given phase distribution, in which case q(x) reduces to a two-dimensional delta function. The magnitude-from-phase problem is thus “less unique” than the phase-from-magnitude problem. However, when the minimum support constraint is known to apply, algorithms of the Fienup-type[4] recover I F(u) I from phase {F(u)} much more effectively than vice versa, and they can handle complexvalued images as efficiently as positive images (which is not true for the phase-from-magnitude problem). The phase must be over-sampled by comparison with what is necessary when both I F(u) I and phase {F(u)} are given, just as the magnitude must be oversampled for the phase-from-magnitude problem, in order that the reconstructed image can be forced to zero outside the minimum support. We feel that the concept of the zero-sheet provides an illuminating geometrical interpretation of the magniproblem, tude(phase)-from-phase( magnitude) which nicely complements Hayes’s detailed algebraic arguments [ 21. From a pragmatic algorithmic viewpoint, the main contribution of this paper is the strategy (described in detail below) for recovering the Fourier magnitude from the modulo II Fourier phase. When the image support is unknown a priori, but the minimum support constraint is known to apply, the Fienup-type algorithms are applied with successively decreasing image-boxes (an image-box is that rectangular region of image-space, centred at the origin with sides parallel to arbitrarily chosen Cartesian axes, outside of which the image is set to zero) until a large increase in the pixel level outside the support is observed. This indicates that the support used previously is the smallest support consistent with the data. Given w(u), as defined by (4), rather than phase{F(u)} itself, we have found the following approach to reconstucting IF(u) I to be effective. The iterations are started by computing a pseudo-random positive image e(x) which vanishes outside a chosen image-box b(x). We next compute the visibility of e(r), which is conveniently written as A(u)exp[i@u)J, where A(u) is positive and O(u) is real. O(u) is then subtracted (modulo2lr) from w(u). At each pixel for which I v(u) -13(u) I is less or
Volume 63, number
1
Fig. 1. Simple image consisting
OPTICS
COMMUNICATIONS
of five isolated pixels.
greater, respectively, than n/2 we redefine 19(u) as w(u) or (w(u) + n) . After recomputing an updated e(x) by Fourier transformation, we set to zero its negative pixels inside, and all of its pixels outside, b(x). An iterative loop is set up thereby. We assess the rate of convergence by calculating, at each iteration, the fraction e of ]e(r) I for which e(x) lies outside b(x) and is negative inside b(r). This procedure is repeated for image-boxes of decreasing size, until the condition noted in the previous paragraph is observed. Fig. 1 shows a simple image consisting of five isolated pixels. Figs. 2 and 3 show the images reconstructed when the linear dimensions of b(x) are set 2OWlarger and smaller, respectively, than those of the correct image-box, which is the b(x) just enclosing the f(r) shown in fig. 1. Fig. 4 shows the variation of 6 with N (the number of iterations) for these image-boxes and the correct image-box. Fig. 2 shows the blurring of f(x) which is inevitable when the image-box is larger than the correct b(x) - note that each pixel is similarly blurred, i.e. it is convolved with a q(x) as defined above. The curve of E versus N for the largest b(x) confirms that convergence is achieved, just as satisfactorily as for the correct b(x). We note that the reconstruction is effectively identical to fig. 1 when b(x) equals the correct image-box. When b(x) is smaller than the correct b(x), fig. 4 indicates that convergence is never achieved, as is also suggested by the sporadic appearance of fig. 3.
1 July 1987
Fig, 2. Reconstruction of image shown in fig. 1 when the linear dimensions of b(x) are 20% larger than those of correct b(x).
These results, together with other computational experience persuade us that we can estimate the linear dimensions of the correct b(x) from inspection of the curves of t versus N. Fig. 5 shows an image containing intricate detail. We added to the oversampled Fourier transform uncorrelated pseudo-random noise, uniformly distributed in the range [ 0, 2n] in phase and [ 0, m] in magnitude where m was 5Okof the DC level. The noise E(u) was forced to be conjugate-symmetric, implying that E( - U) =E* ( U) , thereby ensuring that
Fig. 3. Reconstruction of image shown in fig. 1 when the linear dimensions of b(x) are 20% smaller than those of correct b(x).
13
Volume 63, number
1
OPTICS
the image noise was real. The magnitude of the Fourier transform was discarded and its modulo II phase was retained. Fig. 6 shows the image reconstructed after 500 iterations (when t has decreased to 1.1 x 1Oe3) when b(x) was made equal to the correct image-box. Given the significant noise level and the consequent degradation of high spatial frequency information, fig. 6 suggests that Fienup-type algorithms for recovering magnitude from phase (modulo n) are practical and robust in the presence of noise. This result is comparable to those obtained
Fig. 6. Reconstruction (sampled twice as finely as demanded by the sampling theorem in both directions in Fourier space) to which noise had been added as prescribed in final paragraph of text.
with the conventional Fienup algorithms when applied to the recovery of Fourier phase distributions from noisy magnitude data [ 8,9].
References
Fig. 5.64 x 64 pixel positive image, quantised levels from black to white.
14
to 256 equi-spaced
[ I] A.V. Oppenheim and J.S. Lim, Proc. IEEE 69 (198 1) 529. [2] M.H.Hayes, IEEETrans. ASSP-30 (1982)140. [3] R.H.T. Bates and M.J. McDonnell, Image restoration and reconstruction (Clarendon Press, Oxford, 1986). [4] J.R. Fienup, Appl. Optics 21 (1982) 2758. [S] R.G. Lane, W.R. Fright and R.H.T. Bates, Direct phase retrieval, to appear in IEEE Trans. ASSP-35 (April 1987). [ 61 R.G. Lane and R.H.T. Bates, Automatic multi-dimensional deconvolution, to appear in J.Opt. Sot. Am. A. 4 (January (987). [7] R.H.T. Bates and W.R. Fright, J. Opt. Sot. Am. 73 (1983) 358. [ 81 R.H.T. Bates and D. Mnyama, in: Advances in electronics and electron physics, Vol. 67, ed. P.W. Hawkes (Academic Press, 1986) p. 1. [9] J.R. Fienup, in: Indirect imaging, ed. J.A. Roberts (Cambridge University Press, 1984) p. 99.