- Email: [email protected]
Practical issues on invariant image averaging using the Bispectrum
SIGNAL
PROCESSING
.llp"~a~'.)~.' nlU ELSEVIER
Signal Processing 40 (1994) 119- 128
Practical issues on invariant image averaging using the Bispectrum R. Marabini, J.M. Carazo* Centro Nacional de Biotecnolog~a ( CSIC ). Universidad Autonoma, Canto Blanco, 28049 Madrid, Spain
Received 23 October 1992; revised 6 August 1993 and 5 May 1994
Abstract In this work we present experimental results on the use of the Bispectrum for translational invariant image averaging. The main issue is to systematically explore the applicability of Bispectral averaging to practical cases where the starting images are rather noisy. The main conclusion is that this type of averaging process presents a number of intrinsic instabilities that makes the whole approach difficult to use. We start by proving that the sum of two Bispectra is not necessarily a Bispectrum. We then present a number of practical studies on how many Bispectra from noisy realizations of a given image are needed for a meaningful image recovery as a function of the image signal-to-noise ratio. The question of averaging a non-homogeneous data set is also addressed. Some practical instabilities associated with recursive recovery methods are discussed, as well as some practical clues to overcome them.
Zusammenfassung In dieser Arbeit stellen wir experimentelle Ergebnisse fiber die Verwendung des Bispektrums ffir verschiebungsinvariante Bildmittelung vor. Die Hauptaufgabe besteht darin, die Anwendbarkeit der bispektralen Mittelung in praktischen F/illen systematisch zu untersuchen, wo die Anfangsbilder ziemlich verrauscht sind. Die Hauptschlul3folgerung ist, dab diese Art des Mittelungsvorganges eine Anzahl von intrinsischen Instabilit/iten aufweist, die bewirken, dab die gesamte Methode schwierig zu verwenden ist. Wir beginnen mit dem Beweis, dab die Summe zweier Bispektren nicht notwendigerweise ein Bispektrum ist. Dann stellen wir eine Anzahl von praktischen Untersuchungen darfiber vor, wieviele Bispektren von verrauschten Realisierungen eines gegebenen Bildes ffir eine sinnvolle Bildwiederherstellung als Funktion des Signal-zu-Rauschleistungsverhfiltnisses ben6tigt werden. Die Frage der Mittelung eines nicht-homogenen Datensatzes wird ebenfalls behandelt. Einige praktische Instabilit/iten im Zusammenhang mit rekursiven Wiederherstellungsmethoden werden diskutiert, ebenso wie einige praktische Tricks, um sie zu fiberwinden. Resume Dans cet article, nous presentons des r~sultats exp~rimentaux concernant l'utilisation du bispectre pour le moyennage d'image invariant par translation. Le but principal de ce travail est d'explorer systfmatiquement les possibilitfs d'application du moyennage bispectral clans des cas pratiques, off les images de ddpart sont passablement bruitfes. La conclusion principale est que ce type de moyennage pr+sente des instabilitfs intrinsfques, qui rendent cette approche difficile fi utiliser. Nous commen~ons par dfmontrer que la somme de deux bispectres n'est pas ndcessairement un bispectre. Plusieurs +tudes pratiques sont ensuite pr~sentfes, qui montrent combien de bispectres sont nfcessaires pour
* Corresponding author• Tel.: + 34 1 585 45 43, + 34 1 585 45 10. Fax: + 34 1 5854506. e-mail: [email protected], Carazo @samba.cnb.uam.es. 0165-1684/94/$7.00 (C) 1994 Elsevier Science B.V. All rights reserved SSDl 0 1 6 5 - 1 6 8 4 ( 9 4 ) 0 0 0 6 8 - B
120
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
retrouver une image signifiante ~ipartir de rralisations bruitres de cette image, ceci sous forme d'une fonction du rapport signal-sur-bruit de l'image. On s'intrresse 6galement h la question du moyennage d'un ensemble de donnres nonhomogrnes. Quelques instabilitrs pratiques associres aux mrthodes de recouvrement rrcursives sont discutres, de mrme quelques indications pratiques pour en venir ~ bout sont donnbes. Keywords: Image avering; Bispectrum analysis; Noise issues
1. Introduction There are a number of important practical situations where the final goal is to obtain a good estimation of an a priori unknown object but all the starting data we have is a collection of very noisy frames where the original object is known to be fully contained, although arbitrarily shifted from frame to frame. Issues like image transmission through a noisy channel [I], reduction of coherent speckle noise [5], single-particle electron microscopy of biological specimens [2] or some forms of target identification are clear examples of this problem. Obvious image processing approaches used to enhance the signal-to-noise ratio (SNR) of the final estimate, such as averaging over different frames, are not possible since the object is shifted from image to image. Since we do not precisely know the shape of the object, we cannot design a matched filter approach that would enable us to first of all detect the object in the different images, shift the images, and then average them. Averaging in the Bispectrum domain has recently been proposed as a way to solve this situation [1, 5 7]. The idea is based on the shift-invariance property of the Bispectrum. Starting from a set of noisy images, the Bispectrum of each image is calculated, then averaged, and an improved image is hopefully obtained from the averaged Bispectrum. Throughout this work we assume that images are corrupted by additive white Gaussian noise and that the whole unknown object is always fully contained within the frame. That is, if the original image is noted asf(x), the additive Gaussian noise as n (x), and the corrupted ('experimental') image as g (x), then g(x) = f ( x ) + n(x).
We are going to follow a rather pragmatic and application-oriented point of view. In this way we will address a number of practical problems that arise when the Bispectra of a set of noisy images are to be averaged, such as some of the intrinsic instabilities associated with the Bispectra averaging process itself, as well as determining the minimum number of images that have to be Bispectrum-averaged before a good estimation of the original signal is recovered. The work is organized in the following way. In Section 2 a brief overview of Bispectral analysis is presented. Section 3 is devoted to Bispectral averaging issues, discussing the asymptotic behavior of this process as well as some intrinsic instabilities of the approach in real (finite and discrete) cases. In Section 4 we present practical results on noisy computer-generated images as a function of the image signal-to-noise ratio. Finally, conclusions are presented in Section 5.
2. Overview of Bispectral analysis Denoting a given image by f ( x ) and its corresponding Fourier transform by F(u), then the Bispectrum off(x), F3(ul, u2), is the third moment of F(u) defined as F3(Ul, u2) = F(Ul )F(u2)F( - ul, - u2).
(1)
It can be proven that F3(Ul, u2) is the Fourier transform of the triple autocorrelation off(u). The main properties of the Bispectrum relevant to this work are: (1) The Bispectrum retains information on both magnitude and phase of the original imagef(x). It should be noted that this is a net difference with respect to second-order moments, such as the power spectrum, that lose all phase information.
121
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119 128
(2) The Bispectrum is shift-invariant. That is, the Bispectra of f ( x ) and f ( x + x o ) are the same. (3) The Bispectrum is robust with respect to additive Gaussian noise contamination, meaning that the Bispectrum of a Gaussian process is zero. (4) The Bispectrum is invertible. That is, it is possible to uniquely 'reconstruct'f (x) (up to a linear shift) from F3(ul, u2).
form
B3. . . . .
ged(Ul,
1 " U2) = ~ i_~1 G3,i (Ul, t/2)
= ~ [ F ( u l ) F ( u 2 ) F ( - ut, - u2) i
1
+ F(ul)F(u2)Ni( - ul, - u2) + F(ul)Ni(u2)F( - Ul,
I/2)
- -
+ Ni(ul)F(u2)F( - ul, - u2) + F(ul)Ni(u2)Ni(-Ul,
3. On Bispeetral averaging
+ Ni(ua)F(uz)Ni( - ul, - u2)
The approach that we use in this work is directly related to the Bispectrum properties that we have just mentioned. First of all, the fact that the Bispectrum does not lose the signalphase information and, furthermore, is invertible, makes possible to process the signal in any of the two spaces in a reversible way. Also, the shift-invariance of the Bispectrum indicates a possibly interesting application with regard to those problems where an a priori unknown object is randomly shifted within a frame in a noisy environment and the identification of the object is the main goal. Central to the Bispectrum averaging approach used in this work is the fact that if both the original image f ( x ) and noise n(x) are of zero mean, then the average of an infinite number of noisy Bispectra is an unbiased estimate of the Bispectrum of the original image [3, 8], i.e., 1
lim n ~
n
~
G3,i(ul,
U2) = F 3 ( U l , U2).
-u2)
(2)
n i=
That this is indeed the case can be easily shown in the following way: let gi(x) be the ith experimental image, f ( x ) the original image and hi(X) the ith realization of the noise, then gi(x ) = f (x) + hi(X). Their corresponding Fourier transforms are denoted by Gi(u), F(u) and Ni(u), F 3 ( u l , u 2 ) is the Bispectrum of f ( x ) calculated with Eq. (1). Then the average of n such Bispectra, noted a s B 3 . . . . . ged(/,gl, u2), will take the
+ Ni(ul)Ni(uz)F(
- - //1, - -
112}
+ N i ( u l ) N i ( u 2 ) N i ( - Ul, -- U:,)]
= SSS + SSN + SNS + NSS + SNN + NSN + NNS + NNN.
(3)
We have noted each of the eight terms of Eq. (3) by the pnemotecnic SSS, SSN, and so forth, so as to differentiate among them while indicating how many times either the signal or the noise terms enter into the calculation. Clearly, SSS is the Bispectrum of the original image f ( x ) and N N N is the Bispectrum of the Gaussian noise n(x), that is zero. For the other six terms it is clear that they all tend to zero as n tends to infinity and fmean(X) F ( 0 , 0 ) = n . . . . (x) = N ( 0 , 0 ) - - 0 . However, if n is not infinity, i.e. in all practical cases, some instability problems may occur. In the following we are going to prove that given two noisy realizations g l ( X ) and g2(x) of the same original signal f ( x ) , the averaging of their corresponding Bispectra is not, in general, a Bispectrum, i.e., there is no i m a g e f ' ( x ) such that its Bispectrum are the average of F3, l(/g, V) and F3.2 (u, v). =
Theorem. The addition of two Bispectra is not generally a Bispectrum. Proof. Let F3 be the Bispectrum, defined by (1), where Ul = (h, k) and u2 = (1, m): F3(h, k, l, m) = F(h, k)F(l, m ) F * ( h + 1, k + m).
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
122
Let us assume that there are two D F T s that are identical except at one point.1 We will refer to them as A(h, k) and B(h, k). We are going to divide the proof of the theorem in two parts. For the first one we require the terms A(0, 0) and/or B(0, 0) to be different from zero, while we waive this condition in the second part. Let us further assume for this first part that there is at least one point (r/, K, 2, #) at which A3(r/,K, 2, #) :~ B3(r/, K, 2,#) and (~/,t¢), (2,/J):A (0, 0). Note that this last condition holds true in most applications and that it could only be violated if the Bispectrum were very sparse, in which case it might even be impossible to uniquely reconstruct an image from it. We wish to prove that the following identity is false: A3(h, k, 1, m) + B3(h, k, l, m) = C3(h, k, l, m). That is, there is no D F T C(h, k) that satisfies C3(h, k, l, m) = C(h, k)C(l, m)C*(h + l, k + m). Let us begin taking (h, k) = (0, 0). Then h 3 ( 0 , 0, 1, m) = A(0, O)Z(l, m)A*(l, m)
= A(O, O)[A(I, m)[ 2, B3(0, 0, l, m) = B(0, O)B(1, m)B*(l, m) = B(0, 0)[B(I, m)l 2. If the point at which A and B are different is not (0, 0), 2 then A(0, 0) = B(0, 0) = 2-1/3C(0, 0), and it is possible to divide by A(0, 0) obtaining A3 . . . . . lized(0, 0, l, m)
= A(I, m)A*(l, m) = IA(I, m)l:,
1This assertion could be ambiguous in the following sense. Since the DFT of a real function is antisymmetric, we always store and work with half of it. Therefore, the fact that two DFTs are identical except at one point means that the nonredundant part is identical except at one point. :We are not going to consider the case in which the DC terms of A and B are not the same, although extension of the proof is straightforward.
B3 . . . . . lized(0, 0, I, m) = B(I, m)B*(l, m) = IB(l,
m)l g,
C3 . . . . . lized(0, 0, l, m)
= 21/3C(1, m)C*(l, m) = 21/31C(l, m)l 2, In this way we arrive at the general equation
IC(I, m)l 2 = 2-1/3([A(1, m)l 2 + IB(I, m)12).
(4)
Let us now consider the quartet (r/, ~c, 2, p) such that neither 0/, ~c) nor (2, #) are (0, 0) and the point (r/+ 2, ~c + p) is precisely the point in which A and B are different. Then A3(r/, K, ~,,/a) + B3(r/, K, 2 , / J )
= A(rl, K)A(2, p)A*(q + 2, x + #) + B(r/, x)B(2, p)B*(rl + 2, ~ + #) = A(~/, K)A(2,/~)(A*(r/+ 2, K + #)
+ B*(~ + 2, ~ + ~)).
(5)
If there is a C(h, k) such that the following expression were true: C3(h, k, l, m) = A3(h, k, l, m) + B3(h, k, l, m),
(6)
then Eq. (5) should be equal to
C (7, K)C(2, ~)C*(~ + 2, K + ~) = A(n,tc)A(2, p)(A*(~l + 2, x + p) + B*(~/+ 2, x + p))
(7)
taking moduli on both sides of (7) and using (4), we would then have 3 2[IA(~/+ 2, ~ + ~)12 + IB(r/+ 2, • + #)12 ] = IA(~/+ 2, K + #)l 2 IB(~/+ 2, m +/J)l 2 + 2lA(r/+ 2, x + #)1 IB(~/+ 2, K + #)1 × cos(ang(A(r/+ 2, x + p),B(rl + 2, x + #)), which is the same as (A -- B) 2 = 0, implying that A = B, which is false because we started from the assumption that A and B were different at one point.
3The simplification involves the division by A(h, k)A(1, m) which is always possible because we have assumed that it exits a point which is not zero.
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119 128
Now we are going to deal with the case A(0, 0) = B(0, 0) = 0. The philosophy of the proof is the same as above. First we need to establish a relationship between the moduli similar to the one obtained in Eq. (4). This relationship may be obtained, as in the previous part, if the Bispectrum is not very sparse. 4 []
123
that simplifies into C(h, k)C(l, m)C(q, ~c) = A(h, k)A(l, m)(A(q, ~c) + B(q, to)). Implying that for this point too, 1
Ic(,7, ,OI = ~IA(,7, ,~) + B(,7, ,~)1 [] Lemma 1. Let us assume that there are three Bispectra A3, B3 and C3 and that the corresponding Fourier transforms that rendered these Bispectra are denoted by A, B and C. Then, if C3 were the addition of A3 and B3, the following relationship among the moduli of A, B and C holds:
To continue proving the second part of the theorem, let us use the point (q, ~c, N - 2t/, N - 2~c) (where N is the dimension of the D F T and A(q, K) # B(q, K)) C 3 (q, N, N - 2q, N - 2~c)
If(h, k)l = ~ l A ( h ,
V2
k) + B(h, k)l.
(6)
Let us assume that A and B are different at one point (q, x). Then for those values of A 3 and B 3 for which A(r/, ~:) or B(~/, x) do not enter into the calculation it is obvious that A 3 = B 3. Therefore, for those points ( h , k , l , m ) where (h,k), (l, m), (h + l, k + m) # (q, ~:) we have C3(h, k, 1, m) = 2A3(h, k, l, m)
and by direct substitution into the equation defining the Bispectrum from the Fourier Transform of an image (Eq. (1)), we have that IC(h,k)l = 3 ~ l A ( h , k ) + B(h,k)l = 3 ~ l A ( h , k ) l .
= C(q, ~c)C(N - 2q, N - 2•)C*(N - q, N - ~c) = C ( N - 2~1, N - 2K)CZ(r/, K) = A ( N - 2q, N - 2K)(A(r/, x)2 + B(t/, x) 2) taking moduli at both sides: IC(N - 2q, U - 2tc)[lC2(q, ~c)l = IA(N - 2q, N - 2x)[l(A(q, ~c)2 + B(q, x)21 using Lemma 1 we arrive at IA(N - 2q, N - 2~c)1IC2(q, ~c)12-1/3 =
IA(N
-
2tl, N
-
2x)l I(A(t/, ~c)2 + B(t/, ~c)2)[,
that is, ]C2(t/, ~c)12 1/3 = ] (A(t/, ~c)2 + B(q, ~c)21. As IC2(tl, K)] =
IC(tl,
K)] 2 and using Lemma 1 again
IC2(q, ~c)l = IA2(~/, K) + BZ(q, K)I21/3 = [C(r/, K)I 2 For those other points of A3 and B3 for which A(r/, x) or B(r/, x) do play a role by substitution into Eq. (1) (considering the case only in which h + 1 = r/and k + m = x) by substitution into Eq. (1) we have C(h, k )C(l, m)C(tl, ~c) = A(h, k)A(l, m)A(tl, x) + B(h, k)B(l, m)B(t l, x)
4The Bispectrum interrelates each point of the DFT with every point so, in general, it is not very difficult to find this complete system of equations. In the case of very sparse Bispectra the addition of two Bispectra could be a Bispectrum, but it may be impossible to uniquely reconstruct an image from these Bispectra.
= IA(q, x) + B(q, K)Iz2-2/3; in other words, Ial 2 + Inl 2 - 21hi IQI cos(ang(A, B))
= IA - BI 2 = 0. This holds if A = B, which is false because we started from the assumption that A # B. So if two DFTs are identical except at one point, then the addition of their Bispectra is not, in general, a Bispectrum. It is therefore clear that, in general, the sum of a finite number of Bispectra will not be a Bispectrum, meaning that, strictly speaking, an 'inverse' will not exist. However, most of the existing image
124
R. Marabmi, J.M. Carazo / Signal Processing 40 (1994) 119-128
reconstruction from Bispectrum algorithms will always 'reconstruct' an image from those averaged data, necessarily leading to inconsistencies in the results. The question of how serious these inconsistencies will be is not an easy one to tackle. On the one hand it is very much signal and noise dependent, since many terms in Eq. (3) have mixed signal and noise contributions, while on the other hand it also depends heavily on the particular reconstruction method that will be used. In this work we have addressed this problem in a practical way using different sets of computer-generated noisy images.
4. Practical results We first generated an 'original' binary image of 32 x 32 pixels that was formed by six equally spaced wedges (Fig. l(a)). Up to 900 different realizations of a zero mean white Gaussian noise term were then added to a randomly shifted copy of this image. Two main noise conditions were treated. In one case the noise standard deviation was set to 15. The resulting images (Fig. l(b)) present a noisy background but still very clear features, and they will be referred to in the following as the 'low noise case data set'. A much higher standard deviation of a = 40 was used for the second case. The resulting images (Fig. l(c)) will be referred as the 'high-noisecase data set'.
The general strategy used to study the effect of Bispectrum averaging was the same for the lownoise-case data set and for the high-noise one. In each case a total set of 900 images was subdivided into smaller subsets; the Bispectrum of images within each subset was calculated, then averaged, and then an image was 'reconstructed' from the averaged Bispectrum. As far as Bispectrum reconstruction algorithms are concerned, in this work we have mainly used the recursive method of Sadler and Giannakis I-7]. Other algorithms have been studied and rejected because of unwrapping difficulties (most of the recursive methods that only take care of the phase) or extremely demanding use of resources (global least-squares fitting). It should be noted that the method we have chosen in this study avoids any phase unwrapping and is relatively fast. We first tried averaging subsets of ten images each. The results of averaging the first nine subsets are shown in Fig. 2(a), for the low-noise-case data set, and Fig. 2(b) for the high-noise-case data set. It is clear that in the low-noise case the reconstructed images are in all cases very close to the original noise-free image (Fig. l(a)), while results obtained in the high-noise case are of poor quality (Fig. 2(b)). We then proceeded to enlarge the size of each subgroup of high-noise images up to 100 elements each, clearly improving the results over the previous averaging (Fig. 2(c)). The next step was to average and reconstruct subsets in a cumulative
Fig. 1. Test images: (a) binary signal (max = 115, min = - 11, mean = 0); (b) binary signal plus a white Gaussian noise of/~ = 0, cr = 15; (c) binary signal plus a white Gaussian noise of/~ = 0, cr = 40.
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
125
Fig. 2. Results of averaging a noisy image (from a binary signal) in the bispectral domain. The first image of each line is an original unprocessed frame. (a) Bispectral reconstruction from sets of 10 frames each (noise of a = 15); (b) Bispectral reconstruction from sets of 10 frames each (noise of a = 40); (c) Bispectral reconstruction from sets of 100 frames each (noise of ~ = 40), (d) Bispectral reconstruction from sets of 100, 200, etc. frames (noise of o- = 40); (e) Ideal low-pass filtered versions of the same image shown in (d).
way, that is, averaging the first 100 Bispectra, then the second 200 and so forth. The results obtained this way (Fig. 2(d)) were marginally better than the ones previously obtained without adding up the subsets, indicating that about 200 images were already a large enough data set under these highnoise conditions. Still, there were a few occasions where the reconstructed images presented noticeable artefactual differences with respect to the original noise-free image. In most cases these artifacts presented a stripe-like shape, suggesting that some relatively high-frequency components had been wrongly reconstructed by the Sadler and Giannakis algorithm. Indeed this is not a surprising result considering that this algorithm works recursively in Fourier space. In most cases a slight low-pass filter was enough to attenuate these effects. Nevertheless, it was still possible to find some almost 'pathological' cases where the instabilities of the averaging/reconstructing process appeared particularly strong. In this respect it is worth noting that the reconstruction from 900 frames presented very promin-
ent stripes, while the reconstruction from 897 was visually good (Fig. 3). This behavior has also appeared on other data sets with a different number of averaged Bispectra, although in this case it was particularly noticeable. If more Bispectra are added the reconstruction recovers a better appearance. Another aspect that should be considered in practical situations is the stability of the method with respect to violations in the starting assumptions, especially the one requiring that all measured images were noisy realizations of exactly the same original image. This aspect is particularly important in real data applications where either small changes in the recording or the illumination conditions may occur, or even more drastic inconsistencies among the starting data set are possible. To check this point we proceeded to average two noise-free computer generated images, respectively formed by a set of six and seven equally spaced wedges (Figs. 4(a), (b)). Averaging in real space obviously yielded an image in which the two sets of wedges were mixed (Fig. 4(c)). Bispectral averaging, however, rendered an
126
R. Marabini, J.M. Carazo /Signal Processing 40 (1994) 119-128
Fig. 3. Pathological case. Bispectral reconstruction from (a) 897 frames; (b) 898 frames; (c) 899 frames; (d) 900 frames; demonstrating the instability of the process.
Fig. 4. Inversion of an averaged Bispectrum from an heterogeneous data set: (a) and (b) sample images; (c) direct-average image; (d) Bispectrum average plus inversion.
image where this mixing of features from the six-fold and sevenfold starts appeared blurred by noticeable artifacts (Fig. 4(d)). A direct consequence, we believe, of the fact that the average of two Bispectra is not a Bispectrum. Therefore in this case, direct averaging produced far better results than Bispectral averaging. We have repeated all the above-described processes using nonbinary images. In general the resuits are rather worse. Low-pass filtered versions of the binary wedges already used in Figs. 1 and 2 were computed and Gaussian white noises of different standard deviations were added to them (Fig. 5). Their Bispectra were then calculated and successively averaged and the final average was reconstructed using the Sadler and Giannakis algorithm [7]. In general the results are inferior when compared to binary images, as could be expected from the discussion that follows. Additionally, when the sets of images are small (Figs. 5(a) and (b)), the reconstructed images were independent of the noise standard deviation being used, indicating that, for these small data sets of gray-scale images, the noise is the driving process.
In order to understand these results we must consider two points: First, in our computer-generated images, different frames have different noise, but the noise of each frame in the high SNR case set and in the low one is exactly the same (up to a multiplicative factor). Second, the Bispectrum of a binary or quasi binary and zero mean object can be consistently estimated by the Bispectrum of a noisy image that contains the object [7]. That is, if the image is large enough, so that the statistics of the noise are good, then the Bispectrum of a binary noisy image and the Bispectrum of the signal itself tends to be equal. As it has been proven for binary signals, the third moment of the image 5 ~x,~ g(x)g(y)g*(x + y) can be rewritten as the addition of the third moment of
5The Bispectrum is the Fourier transform of the third moment of the image.
R. Marabini. J.M. Carazo /Signal Processing 40 (1994) 119-128
127
Fig. 5. Results of averaging a noisy data set of gray-scale images. The first image of each line is an original unprocessed frame: (a) Bispectral reconstruction from sets of l0 frames each (noise of a = 15);(b) Bispectral reconstruction from sets of 10 frames each (noise of a = 40); (c) Bispectral reconstruction from sets of 100 frames each (noise of a = 40};(d) Bispectral reconstruction from sets of 100, 200, etc. frames (noise of a = 40).
the signal, the third m o m e n t of the noise (that tends to zero as the image gets larger) and terms like
~ x~l,+,ycL,
= J m e a n ( r ( x -- y ) ) =
~
f2 (x)n(x + y)
xcL-,yEL
+
~
f+ (x)f (y)n(x + y)
xEL+,j'eL "2 =J+(x)
Z
n(x+y)
x~l.*,yEL+
+ f2_(x)
~
n(x + y)
.r~L-,y~L
+ Jr+(x).L (y)
~
n(x + y)
x~L+,yeL-
= 0 if
Z
n(x + y) =
xeL+.y~L+
=
~
Z xeL-,y~L-
n(x + y) = 0
x~L+,y~L
and
~f(x)n(y)nlx + y) n(y)n(x + y)
=/+(x) x~l,+,yEL+
n(y)nIx + y) if L+ and L
are large enough
fz+(x)n(x + y)
+
~
xeL+.yeL+
= (f+(x) + f ( x ) ) ( r ( x - y ) )
~ f ( x ) f ( y ) n ( x + y) =
+ f Ix)
n(x + y)
O,
where f+ and J2 are the two possible values of the binary signal, L + and L must be interpreted as the set of points in which the value of the signal is f+ and ~ , respectively, and r ( x - y) is the noise autocorrelation. Obviously, if the signal is not binary but grayscaled we then obtain the same terms for each different pixel value. However, if the n u m b e r of different gray levels is high enough and since we have been working with 32 by 32 pixel images, the number of points that belong to a significant number of Li would be low, meaning that the terms that before tended to zero may now be higher than the signal third moment. Besides, since the signal is the same in both the high-noise case and the low-noise one, and the noise itself is related by a constant factor, as we c o m m e n t e d before, then these terms are now equal in both data sets (except for a multiplicative factor that is automatically compensated by the display algorithm).
128
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
5. Conclusions
In this work we have studied in a rather practical way the feasibility of Bispectral averaging as a way to identify an a priori unknown object from a set of noisy frames presenting the object wholly contained, but otherwise arbitrarily shifted, within the frame. The main conclusions to be derived from this work are as follows. Averaging in Bispectrum domain is an unstable process. The reason being that, in general, the average of a finite number of Bispectra will not be a Bispectrum, i.e. the average will not be invertible. If standard algorithms for image reconstruction from Bispectral data are still used on such averaged data, the resulting image will tend to be closer to the original noise-free image as the number of terms to be averaged increases. The minimum number of images to be averaged depends on the signal-to-noise ratio of the individual frames and on the particular reconstruction algorithm that is being used in a complex and signal-dependent way. However, numbers in the order of the tens of images for low-noise cases and of the hundreds for more serious noise conditions should be expected to be necessary for binary signals. Averaging in Bispectral domain, although unstable, is still able to render an estimation with a high signal-to-noise ratio if the starting data set is composed of binary (or almost binary) images.
However, and within the limitations of this work, it seems that this approach cannot be used in noisy gray-scaled images.
References [1] S.A. Dianat and M.R. Raghuveer, "Fast algorithms for phase and magnitude recovering" Opt. Engrg., Vol. 29, No. 5, May 1990, pp. 504512. [2] J. Frank, A. Verschoor and M. Boublik, "Computer averaging of electron micrographs of 40S ribosomal subunits", Science, Vol. 214, December 1981, pp. 1353-1355. [3] G.B. Giannakis, "Signal reconstruction from multiple correlations: Frequency and time domain approaches", J. Opt. Soc. Amer. A, Vol. 6, No. 5, May 1989, pp. 682~97. [4] C.A. Haniff, "Least-squares Fourier phase estimation from the module 2n Bispectrum phase", J. Opt. Soc. Amer. A., Vol. 8, No. 1, January 1991, pp. 134-140. [5] S. Jim, S. Wear and M.R. Raghuveer, "Reconstruction of speckled images using Bispectra", J. Opt. Soc. Amer. A., Vol. 9, No. 3, March 1992, pp. 371 376. [6] A.W. Lohmann, "Shift and rotation tolerant image recovery from triple correlations", Optik, Vol. 78, No. 3, 1986, pp. 127-131. [7] B.M. Sadler and G. B. Giannakis, "Shift- and rotationinvariant object reconstruction using the Bispectrum", J. Opt. Soc. Amer. A., Vol. 9, No. 1, January 1992, pp. 57q59. [8] G. Sundaramoorthy, M.R. Raghuveer and S.A. Dianat, "Bispectral reconstruction of signals in noise: Amplitude reconstruction issues", IEEE Trans. Acoust. Speech Signal Process., Vol. 38, No. 7, July 1990, pp. 1297-1306. I-9] M.K. Tsatsanis and G.B Giannakis, "Object and texture classification using higher-order statistics", IEEE Trans. Acoust Speech Signal Process., Vol. 38. July 1990, pp. 1284-1296.
PROCESSING
.llp"~a~'.)~.' nlU ELSEVIER
Signal Processing 40 (1994) 119- 128
Practical issues on invariant image averaging using the Bispectrum R. Marabini, J.M. Carazo* Centro Nacional de Biotecnolog~a ( CSIC ). Universidad Autonoma, Canto Blanco, 28049 Madrid, Spain
Received 23 October 1992; revised 6 August 1993 and 5 May 1994
Abstract In this work we present experimental results on the use of the Bispectrum for translational invariant image averaging. The main issue is to systematically explore the applicability of Bispectral averaging to practical cases where the starting images are rather noisy. The main conclusion is that this type of averaging process presents a number of intrinsic instabilities that makes the whole approach difficult to use. We start by proving that the sum of two Bispectra is not necessarily a Bispectrum. We then present a number of practical studies on how many Bispectra from noisy realizations of a given image are needed for a meaningful image recovery as a function of the image signal-to-noise ratio. The question of averaging a non-homogeneous data set is also addressed. Some practical instabilities associated with recursive recovery methods are discussed, as well as some practical clues to overcome them.
Zusammenfassung In dieser Arbeit stellen wir experimentelle Ergebnisse fiber die Verwendung des Bispektrums ffir verschiebungsinvariante Bildmittelung vor. Die Hauptaufgabe besteht darin, die Anwendbarkeit der bispektralen Mittelung in praktischen F/illen systematisch zu untersuchen, wo die Anfangsbilder ziemlich verrauscht sind. Die Hauptschlul3folgerung ist, dab diese Art des Mittelungsvorganges eine Anzahl von intrinsischen Instabilit/iten aufweist, die bewirken, dab die gesamte Methode schwierig zu verwenden ist. Wir beginnen mit dem Beweis, dab die Summe zweier Bispektren nicht notwendigerweise ein Bispektrum ist. Dann stellen wir eine Anzahl von praktischen Untersuchungen darfiber vor, wieviele Bispektren von verrauschten Realisierungen eines gegebenen Bildes ffir eine sinnvolle Bildwiederherstellung als Funktion des Signal-zu-Rauschleistungsverhfiltnisses ben6tigt werden. Die Frage der Mittelung eines nicht-homogenen Datensatzes wird ebenfalls behandelt. Einige praktische Instabilit/iten im Zusammenhang mit rekursiven Wiederherstellungsmethoden werden diskutiert, ebenso wie einige praktische Tricks, um sie zu fiberwinden. Resume Dans cet article, nous presentons des r~sultats exp~rimentaux concernant l'utilisation du bispectre pour le moyennage d'image invariant par translation. Le but principal de ce travail est d'explorer systfmatiquement les possibilitfs d'application du moyennage bispectral clans des cas pratiques, off les images de ddpart sont passablement bruitfes. La conclusion principale est que ce type de moyennage pr+sente des instabilitfs intrinsfques, qui rendent cette approche difficile fi utiliser. Nous commen~ons par dfmontrer que la somme de deux bispectres n'est pas ndcessairement un bispectre. Plusieurs +tudes pratiques sont ensuite pr~sentfes, qui montrent combien de bispectres sont nfcessaires pour
* Corresponding author• Tel.: + 34 1 585 45 43, + 34 1 585 45 10. Fax: + 34 1 5854506. e-mail: [email protected], Carazo @samba.cnb.uam.es. 0165-1684/94/$7.00 (C) 1994 Elsevier Science B.V. All rights reserved SSDl 0 1 6 5 - 1 6 8 4 ( 9 4 ) 0 0 0 6 8 - B
120
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
retrouver une image signifiante ~ipartir de rralisations bruitres de cette image, ceci sous forme d'une fonction du rapport signal-sur-bruit de l'image. On s'intrresse 6galement h la question du moyennage d'un ensemble de donnres nonhomogrnes. Quelques instabilitrs pratiques associres aux mrthodes de recouvrement rrcursives sont discutres, de mrme quelques indications pratiques pour en venir ~ bout sont donnbes. Keywords: Image avering; Bispectrum analysis; Noise issues
1. Introduction There are a number of important practical situations where the final goal is to obtain a good estimation of an a priori unknown object but all the starting data we have is a collection of very noisy frames where the original object is known to be fully contained, although arbitrarily shifted from frame to frame. Issues like image transmission through a noisy channel [I], reduction of coherent speckle noise [5], single-particle electron microscopy of biological specimens [2] or some forms of target identification are clear examples of this problem. Obvious image processing approaches used to enhance the signal-to-noise ratio (SNR) of the final estimate, such as averaging over different frames, are not possible since the object is shifted from image to image. Since we do not precisely know the shape of the object, we cannot design a matched filter approach that would enable us to first of all detect the object in the different images, shift the images, and then average them. Averaging in the Bispectrum domain has recently been proposed as a way to solve this situation [1, 5 7]. The idea is based on the shift-invariance property of the Bispectrum. Starting from a set of noisy images, the Bispectrum of each image is calculated, then averaged, and an improved image is hopefully obtained from the averaged Bispectrum. Throughout this work we assume that images are corrupted by additive white Gaussian noise and that the whole unknown object is always fully contained within the frame. That is, if the original image is noted asf(x), the additive Gaussian noise as n (x), and the corrupted ('experimental') image as g (x), then g(x) = f ( x ) + n(x).
We are going to follow a rather pragmatic and application-oriented point of view. In this way we will address a number of practical problems that arise when the Bispectra of a set of noisy images are to be averaged, such as some of the intrinsic instabilities associated with the Bispectra averaging process itself, as well as determining the minimum number of images that have to be Bispectrum-averaged before a good estimation of the original signal is recovered. The work is organized in the following way. In Section 2 a brief overview of Bispectral analysis is presented. Section 3 is devoted to Bispectral averaging issues, discussing the asymptotic behavior of this process as well as some intrinsic instabilities of the approach in real (finite and discrete) cases. In Section 4 we present practical results on noisy computer-generated images as a function of the image signal-to-noise ratio. Finally, conclusions are presented in Section 5.
2. Overview of Bispectral analysis Denoting a given image by f ( x ) and its corresponding Fourier transform by F(u), then the Bispectrum off(x), F3(ul, u2), is the third moment of F(u) defined as F3(Ul, u2) = F(Ul )F(u2)F( - ul, - u2).
(1)
It can be proven that F3(Ul, u2) is the Fourier transform of the triple autocorrelation off(u). The main properties of the Bispectrum relevant to this work are: (1) The Bispectrum retains information on both magnitude and phase of the original imagef(x). It should be noted that this is a net difference with respect to second-order moments, such as the power spectrum, that lose all phase information.
121
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119 128
(2) The Bispectrum is shift-invariant. That is, the Bispectra of f ( x ) and f ( x + x o ) are the same. (3) The Bispectrum is robust with respect to additive Gaussian noise contamination, meaning that the Bispectrum of a Gaussian process is zero. (4) The Bispectrum is invertible. That is, it is possible to uniquely 'reconstruct'f (x) (up to a linear shift) from F3(ul, u2).
form
B3. . . . .
ged(Ul,
1 " U2) = ~ i_~1 G3,i (Ul, t/2)
= ~ [ F ( u l ) F ( u 2 ) F ( - ut, - u2) i
1
+ F(ul)F(u2)Ni( - ul, - u2) + F(ul)Ni(u2)F( - Ul,
I/2)
- -
+ Ni(ul)F(u2)F( - ul, - u2) + F(ul)Ni(u2)Ni(-Ul,
3. On Bispeetral averaging
+ Ni(ua)F(uz)Ni( - ul, - u2)
The approach that we use in this work is directly related to the Bispectrum properties that we have just mentioned. First of all, the fact that the Bispectrum does not lose the signalphase information and, furthermore, is invertible, makes possible to process the signal in any of the two spaces in a reversible way. Also, the shift-invariance of the Bispectrum indicates a possibly interesting application with regard to those problems where an a priori unknown object is randomly shifted within a frame in a noisy environment and the identification of the object is the main goal. Central to the Bispectrum averaging approach used in this work is the fact that if both the original image f ( x ) and noise n(x) are of zero mean, then the average of an infinite number of noisy Bispectra is an unbiased estimate of the Bispectrum of the original image [3, 8], i.e., 1
lim n ~
n
~
G3,i(ul,
U2) = F 3 ( U l , U2).
-u2)
(2)
n i=
That this is indeed the case can be easily shown in the following way: let gi(x) be the ith experimental image, f ( x ) the original image and hi(X) the ith realization of the noise, then gi(x ) = f (x) + hi(X). Their corresponding Fourier transforms are denoted by Gi(u), F(u) and Ni(u), F 3 ( u l , u 2 ) is the Bispectrum of f ( x ) calculated with Eq. (1). Then the average of n such Bispectra, noted a s B 3 . . . . . ged(/,gl, u2), will take the
+ Ni(ul)Ni(uz)F(
- - //1, - -
112}
+ N i ( u l ) N i ( u 2 ) N i ( - Ul, -- U:,)]
= SSS + SSN + SNS + NSS + SNN + NSN + NNS + NNN.
(3)
We have noted each of the eight terms of Eq. (3) by the pnemotecnic SSS, SSN, and so forth, so as to differentiate among them while indicating how many times either the signal or the noise terms enter into the calculation. Clearly, SSS is the Bispectrum of the original image f ( x ) and N N N is the Bispectrum of the Gaussian noise n(x), that is zero. For the other six terms it is clear that they all tend to zero as n tends to infinity and fmean(X) F ( 0 , 0 ) = n . . . . (x) = N ( 0 , 0 ) - - 0 . However, if n is not infinity, i.e. in all practical cases, some instability problems may occur. In the following we are going to prove that given two noisy realizations g l ( X ) and g2(x) of the same original signal f ( x ) , the averaging of their corresponding Bispectra is not, in general, a Bispectrum, i.e., there is no i m a g e f ' ( x ) such that its Bispectrum are the average of F3, l(/g, V) and F3.2 (u, v). =
Theorem. The addition of two Bispectra is not generally a Bispectrum. Proof. Let F3 be the Bispectrum, defined by (1), where Ul = (h, k) and u2 = (1, m): F3(h, k, l, m) = F(h, k)F(l, m ) F * ( h + 1, k + m).
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
122
Let us assume that there are two D F T s that are identical except at one point.1 We will refer to them as A(h, k) and B(h, k). We are going to divide the proof of the theorem in two parts. For the first one we require the terms A(0, 0) and/or B(0, 0) to be different from zero, while we waive this condition in the second part. Let us further assume for this first part that there is at least one point (r/, K, 2, #) at which A3(r/,K, 2, #) :~ B3(r/, K, 2,#) and (~/,t¢), (2,/J):A (0, 0). Note that this last condition holds true in most applications and that it could only be violated if the Bispectrum were very sparse, in which case it might even be impossible to uniquely reconstruct an image from it. We wish to prove that the following identity is false: A3(h, k, 1, m) + B3(h, k, l, m) = C3(h, k, l, m). That is, there is no D F T C(h, k) that satisfies C3(h, k, l, m) = C(h, k)C(l, m)C*(h + l, k + m). Let us begin taking (h, k) = (0, 0). Then h 3 ( 0 , 0, 1, m) = A(0, O)Z(l, m)A*(l, m)
= A(O, O)[A(I, m)[ 2, B3(0, 0, l, m) = B(0, O)B(1, m)B*(l, m) = B(0, 0)[B(I, m)l 2. If the point at which A and B are different is not (0, 0), 2 then A(0, 0) = B(0, 0) = 2-1/3C(0, 0), and it is possible to divide by A(0, 0) obtaining A3 . . . . . lized(0, 0, l, m)
= A(I, m)A*(l, m) = IA(I, m)l:,
1This assertion could be ambiguous in the following sense. Since the DFT of a real function is antisymmetric, we always store and work with half of it. Therefore, the fact that two DFTs are identical except at one point means that the nonredundant part is identical except at one point. :We are not going to consider the case in which the DC terms of A and B are not the same, although extension of the proof is straightforward.
B3 . . . . . lized(0, 0, I, m) = B(I, m)B*(l, m) = IB(l,
m)l g,
C3 . . . . . lized(0, 0, l, m)
= 21/3C(1, m)C*(l, m) = 21/31C(l, m)l 2, In this way we arrive at the general equation
IC(I, m)l 2 = 2-1/3([A(1, m)l 2 + IB(I, m)12).
(4)
Let us now consider the quartet (r/, ~c, 2, p) such that neither 0/, ~c) nor (2, #) are (0, 0) and the point (r/+ 2, ~c + p) is precisely the point in which A and B are different. Then A3(r/, K, ~,,/a) + B3(r/, K, 2 , / J )
= A(rl, K)A(2, p)A*(q + 2, x + #) + B(r/, x)B(2, p)B*(rl + 2, ~ + #) = A(~/, K)A(2,/~)(A*(r/+ 2, K + #)
+ B*(~ + 2, ~ + ~)).
(5)
If there is a C(h, k) such that the following expression were true: C3(h, k, l, m) = A3(h, k, l, m) + B3(h, k, l, m),
(6)
then Eq. (5) should be equal to
C (7, K)C(2, ~)C*(~ + 2, K + ~) = A(n,tc)A(2, p)(A*(~l + 2, x + p) + B*(~/+ 2, x + p))
(7)
taking moduli on both sides of (7) and using (4), we would then have 3 2[IA(~/+ 2, ~ + ~)12 + IB(r/+ 2, • + #)12 ] = IA(~/+ 2, K + #)l 2 IB(~/+ 2, m +/J)l 2 + 2lA(r/+ 2, x + #)1 IB(~/+ 2, K + #)1 × cos(ang(A(r/+ 2, x + p),B(rl + 2, x + #)), which is the same as (A -- B) 2 = 0, implying that A = B, which is false because we started from the assumption that A and B were different at one point.
3The simplification involves the division by A(h, k)A(1, m) which is always possible because we have assumed that it exits a point which is not zero.
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119 128
Now we are going to deal with the case A(0, 0) = B(0, 0) = 0. The philosophy of the proof is the same as above. First we need to establish a relationship between the moduli similar to the one obtained in Eq. (4). This relationship may be obtained, as in the previous part, if the Bispectrum is not very sparse. 4 []
123
that simplifies into C(h, k)C(l, m)C(q, ~c) = A(h, k)A(l, m)(A(q, ~c) + B(q, to)). Implying that for this point too, 1
Ic(,7, ,OI = ~IA(,7, ,~) + B(,7, ,~)1 [] Lemma 1. Let us assume that there are three Bispectra A3, B3 and C3 and that the corresponding Fourier transforms that rendered these Bispectra are denoted by A, B and C. Then, if C3 were the addition of A3 and B3, the following relationship among the moduli of A, B and C holds:
To continue proving the second part of the theorem, let us use the point (q, ~c, N - 2t/, N - 2~c) (where N is the dimension of the D F T and A(q, K) # B(q, K)) C 3 (q, N, N - 2q, N - 2~c)
If(h, k)l = ~ l A ( h ,
V2
k) + B(h, k)l.
(6)
Let us assume that A and B are different at one point (q, x). Then for those values of A 3 and B 3 for which A(r/, ~:) or B(~/, x) do not enter into the calculation it is obvious that A 3 = B 3. Therefore, for those points ( h , k , l , m ) where (h,k), (l, m), (h + l, k + m) # (q, ~:) we have C3(h, k, 1, m) = 2A3(h, k, l, m)
and by direct substitution into the equation defining the Bispectrum from the Fourier Transform of an image (Eq. (1)), we have that IC(h,k)l = 3 ~ l A ( h , k ) + B(h,k)l = 3 ~ l A ( h , k ) l .
= C(q, ~c)C(N - 2q, N - 2•)C*(N - q, N - ~c) = C ( N - 2~1, N - 2K)CZ(r/, K) = A ( N - 2q, N - 2K)(A(r/, x)2 + B(t/, x) 2) taking moduli at both sides: IC(N - 2q, U - 2tc)[lC2(q, ~c)l = IA(N - 2q, N - 2x)[l(A(q, ~c)2 + B(q, x)21 using Lemma 1 we arrive at IA(N - 2q, N - 2~c)1IC2(q, ~c)12-1/3 =
IA(N
-
2tl, N
-
2x)l I(A(t/, ~c)2 + B(t/, ~c)2)[,
that is, ]C2(t/, ~c)12 1/3 = ] (A(t/, ~c)2 + B(q, ~c)21. As IC2(tl, K)] =
IC(tl,
K)] 2 and using Lemma 1 again
IC2(q, ~c)l = IA2(~/, K) + BZ(q, K)I21/3 = [C(r/, K)I 2 For those other points of A3 and B3 for which A(r/, x) or B(r/, x) do play a role by substitution into Eq. (1) (considering the case only in which h + 1 = r/and k + m = x) by substitution into Eq. (1) we have C(h, k )C(l, m)C(tl, ~c) = A(h, k)A(l, m)A(tl, x) + B(h, k)B(l, m)B(t l, x)
4The Bispectrum interrelates each point of the DFT with every point so, in general, it is not very difficult to find this complete system of equations. In the case of very sparse Bispectra the addition of two Bispectra could be a Bispectrum, but it may be impossible to uniquely reconstruct an image from these Bispectra.
= IA(q, x) + B(q, K)Iz2-2/3; in other words, Ial 2 + Inl 2 - 21hi IQI cos(ang(A, B))
= IA - BI 2 = 0. This holds if A = B, which is false because we started from the assumption that A # B. So if two DFTs are identical except at one point, then the addition of their Bispectra is not, in general, a Bispectrum. It is therefore clear that, in general, the sum of a finite number of Bispectra will not be a Bispectrum, meaning that, strictly speaking, an 'inverse' will not exist. However, most of the existing image
124
R. Marabmi, J.M. Carazo / Signal Processing 40 (1994) 119-128
reconstruction from Bispectrum algorithms will always 'reconstruct' an image from those averaged data, necessarily leading to inconsistencies in the results. The question of how serious these inconsistencies will be is not an easy one to tackle. On the one hand it is very much signal and noise dependent, since many terms in Eq. (3) have mixed signal and noise contributions, while on the other hand it also depends heavily on the particular reconstruction method that will be used. In this work we have addressed this problem in a practical way using different sets of computer-generated noisy images.
4. Practical results We first generated an 'original' binary image of 32 x 32 pixels that was formed by six equally spaced wedges (Fig. l(a)). Up to 900 different realizations of a zero mean white Gaussian noise term were then added to a randomly shifted copy of this image. Two main noise conditions were treated. In one case the noise standard deviation was set to 15. The resulting images (Fig. l(b)) present a noisy background but still very clear features, and they will be referred to in the following as the 'low noise case data set'. A much higher standard deviation of a = 40 was used for the second case. The resulting images (Fig. l(c)) will be referred as the 'high-noisecase data set'.
The general strategy used to study the effect of Bispectrum averaging was the same for the lownoise-case data set and for the high-noise one. In each case a total set of 900 images was subdivided into smaller subsets; the Bispectrum of images within each subset was calculated, then averaged, and then an image was 'reconstructed' from the averaged Bispectrum. As far as Bispectrum reconstruction algorithms are concerned, in this work we have mainly used the recursive method of Sadler and Giannakis I-7]. Other algorithms have been studied and rejected because of unwrapping difficulties (most of the recursive methods that only take care of the phase) or extremely demanding use of resources (global least-squares fitting). It should be noted that the method we have chosen in this study avoids any phase unwrapping and is relatively fast. We first tried averaging subsets of ten images each. The results of averaging the first nine subsets are shown in Fig. 2(a), for the low-noise-case data set, and Fig. 2(b) for the high-noise-case data set. It is clear that in the low-noise case the reconstructed images are in all cases very close to the original noise-free image (Fig. l(a)), while results obtained in the high-noise case are of poor quality (Fig. 2(b)). We then proceeded to enlarge the size of each subgroup of high-noise images up to 100 elements each, clearly improving the results over the previous averaging (Fig. 2(c)). The next step was to average and reconstruct subsets in a cumulative
Fig. 1. Test images: (a) binary signal (max = 115, min = - 11, mean = 0); (b) binary signal plus a white Gaussian noise of/~ = 0, cr = 15; (c) binary signal plus a white Gaussian noise of/~ = 0, cr = 40.
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
125
Fig. 2. Results of averaging a noisy image (from a binary signal) in the bispectral domain. The first image of each line is an original unprocessed frame. (a) Bispectral reconstruction from sets of 10 frames each (noise of a = 15); (b) Bispectral reconstruction from sets of 10 frames each (noise of a = 40); (c) Bispectral reconstruction from sets of 100 frames each (noise of ~ = 40), (d) Bispectral reconstruction from sets of 100, 200, etc. frames (noise of o- = 40); (e) Ideal low-pass filtered versions of the same image shown in (d).
way, that is, averaging the first 100 Bispectra, then the second 200 and so forth. The results obtained this way (Fig. 2(d)) were marginally better than the ones previously obtained without adding up the subsets, indicating that about 200 images were already a large enough data set under these highnoise conditions. Still, there were a few occasions where the reconstructed images presented noticeable artefactual differences with respect to the original noise-free image. In most cases these artifacts presented a stripe-like shape, suggesting that some relatively high-frequency components had been wrongly reconstructed by the Sadler and Giannakis algorithm. Indeed this is not a surprising result considering that this algorithm works recursively in Fourier space. In most cases a slight low-pass filter was enough to attenuate these effects. Nevertheless, it was still possible to find some almost 'pathological' cases where the instabilities of the averaging/reconstructing process appeared particularly strong. In this respect it is worth noting that the reconstruction from 900 frames presented very promin-
ent stripes, while the reconstruction from 897 was visually good (Fig. 3). This behavior has also appeared on other data sets with a different number of averaged Bispectra, although in this case it was particularly noticeable. If more Bispectra are added the reconstruction recovers a better appearance. Another aspect that should be considered in practical situations is the stability of the method with respect to violations in the starting assumptions, especially the one requiring that all measured images were noisy realizations of exactly the same original image. This aspect is particularly important in real data applications where either small changes in the recording or the illumination conditions may occur, or even more drastic inconsistencies among the starting data set are possible. To check this point we proceeded to average two noise-free computer generated images, respectively formed by a set of six and seven equally spaced wedges (Figs. 4(a), (b)). Averaging in real space obviously yielded an image in which the two sets of wedges were mixed (Fig. 4(c)). Bispectral averaging, however, rendered an
126
R. Marabini, J.M. Carazo /Signal Processing 40 (1994) 119-128
Fig. 3. Pathological case. Bispectral reconstruction from (a) 897 frames; (b) 898 frames; (c) 899 frames; (d) 900 frames; demonstrating the instability of the process.
Fig. 4. Inversion of an averaged Bispectrum from an heterogeneous data set: (a) and (b) sample images; (c) direct-average image; (d) Bispectrum average plus inversion.
image where this mixing of features from the six-fold and sevenfold starts appeared blurred by noticeable artifacts (Fig. 4(d)). A direct consequence, we believe, of the fact that the average of two Bispectra is not a Bispectrum. Therefore in this case, direct averaging produced far better results than Bispectral averaging. We have repeated all the above-described processes using nonbinary images. In general the resuits are rather worse. Low-pass filtered versions of the binary wedges already used in Figs. 1 and 2 were computed and Gaussian white noises of different standard deviations were added to them (Fig. 5). Their Bispectra were then calculated and successively averaged and the final average was reconstructed using the Sadler and Giannakis algorithm [7]. In general the results are inferior when compared to binary images, as could be expected from the discussion that follows. Additionally, when the sets of images are small (Figs. 5(a) and (b)), the reconstructed images were independent of the noise standard deviation being used, indicating that, for these small data sets of gray-scale images, the noise is the driving process.
In order to understand these results we must consider two points: First, in our computer-generated images, different frames have different noise, but the noise of each frame in the high SNR case set and in the low one is exactly the same (up to a multiplicative factor). Second, the Bispectrum of a binary or quasi binary and zero mean object can be consistently estimated by the Bispectrum of a noisy image that contains the object [7]. That is, if the image is large enough, so that the statistics of the noise are good, then the Bispectrum of a binary noisy image and the Bispectrum of the signal itself tends to be equal. As it has been proven for binary signals, the third moment of the image 5 ~x,~ g(x)g(y)g*(x + y) can be rewritten as the addition of the third moment of
5The Bispectrum is the Fourier transform of the third moment of the image.
R. Marabini. J.M. Carazo /Signal Processing 40 (1994) 119-128
127
Fig. 5. Results of averaging a noisy data set of gray-scale images. The first image of each line is an original unprocessed frame: (a) Bispectral reconstruction from sets of l0 frames each (noise of a = 15);(b) Bispectral reconstruction from sets of 10 frames each (noise of a = 40); (c) Bispectral reconstruction from sets of 100 frames each (noise of a = 40};(d) Bispectral reconstruction from sets of 100, 200, etc. frames (noise of a = 40).
the signal, the third m o m e n t of the noise (that tends to zero as the image gets larger) and terms like
~ x~l,+,ycL,
= J m e a n ( r ( x -- y ) ) =
~
f2 (x)n(x + y)
xcL-,yEL
+
~
f+ (x)f (y)n(x + y)
xEL+,j'eL "2 =J+(x)
Z
n(x+y)
x~l.*,yEL+
+ f2_(x)
~
n(x + y)
.r~L-,y~L
+ Jr+(x).L (y)
~
n(x + y)
x~L+,yeL-
= 0 if
Z
n(x + y) =
xeL+.y~L+
=
~
Z xeL-,y~L-
n(x + y) = 0
x~L+,y~L
and
~f(x)n(y)nlx + y) n(y)n(x + y)
=/+(x) x~l,+,yEL+
n(y)nIx + y) if L+ and L
are large enough
fz+(x)n(x + y)
+
~
xeL+.yeL+
= (f+(x) + f ( x ) ) ( r ( x - y ) )
~ f ( x ) f ( y ) n ( x + y) =
+ f Ix)
n(x + y)
O,
where f+ and J2 are the two possible values of the binary signal, L + and L must be interpreted as the set of points in which the value of the signal is f+ and ~ , respectively, and r ( x - y) is the noise autocorrelation. Obviously, if the signal is not binary but grayscaled we then obtain the same terms for each different pixel value. However, if the n u m b e r of different gray levels is high enough and since we have been working with 32 by 32 pixel images, the number of points that belong to a significant number of Li would be low, meaning that the terms that before tended to zero may now be higher than the signal third moment. Besides, since the signal is the same in both the high-noise case and the low-noise one, and the noise itself is related by a constant factor, as we c o m m e n t e d before, then these terms are now equal in both data sets (except for a multiplicative factor that is automatically compensated by the display algorithm).
128
R. Marabini, J.M. Carazo / Signal Processing 40 (1994) 119-128
5. Conclusions
In this work we have studied in a rather practical way the feasibility of Bispectral averaging as a way to identify an a priori unknown object from a set of noisy frames presenting the object wholly contained, but otherwise arbitrarily shifted, within the frame. The main conclusions to be derived from this work are as follows. Averaging in Bispectrum domain is an unstable process. The reason being that, in general, the average of a finite number of Bispectra will not be a Bispectrum, i.e. the average will not be invertible. If standard algorithms for image reconstruction from Bispectral data are still used on such averaged data, the resulting image will tend to be closer to the original noise-free image as the number of terms to be averaged increases. The minimum number of images to be averaged depends on the signal-to-noise ratio of the individual frames and on the particular reconstruction algorithm that is being used in a complex and signal-dependent way. However, numbers in the order of the tens of images for low-noise cases and of the hundreds for more serious noise conditions should be expected to be necessary for binary signals. Averaging in Bispectral domain, although unstable, is still able to render an estimation with a high signal-to-noise ratio if the starting data set is composed of binary (or almost binary) images.
However, and within the limitations of this work, it seems that this approach cannot be used in noisy gray-scaled images.
References [1] S.A. Dianat and M.R. Raghuveer, "Fast algorithms for phase and magnitude recovering" Opt. Engrg., Vol. 29, No. 5, May 1990, pp. 504512. [2] J. Frank, A. Verschoor and M. Boublik, "Computer averaging of electron micrographs of 40S ribosomal subunits", Science, Vol. 214, December 1981, pp. 1353-1355. [3] G.B. Giannakis, "Signal reconstruction from multiple correlations: Frequency and time domain approaches", J. Opt. Soc. Amer. A, Vol. 6, No. 5, May 1989, pp. 682~97. [4] C.A. Haniff, "Least-squares Fourier phase estimation from the module 2n Bispectrum phase", J. Opt. Soc. Amer. A., Vol. 8, No. 1, January 1991, pp. 134-140. [5] S. Jim, S. Wear and M.R. Raghuveer, "Reconstruction of speckled images using Bispectra", J. Opt. Soc. Amer. A., Vol. 9, No. 3, March 1992, pp. 371 376. [6] A.W. Lohmann, "Shift and rotation tolerant image recovery from triple correlations", Optik, Vol. 78, No. 3, 1986, pp. 127-131. [7] B.M. Sadler and G. B. Giannakis, "Shift- and rotationinvariant object reconstruction using the Bispectrum", J. Opt. Soc. Amer. A., Vol. 9, No. 1, January 1992, pp. 57q59. [8] G. Sundaramoorthy, M.R. Raghuveer and S.A. Dianat, "Bispectral reconstruction of signals in noise: Amplitude reconstruction issues", IEEE Trans. Acoust. Speech Signal Process., Vol. 38, No. 7, July 1990, pp. 1297-1306. I-9] M.K. Tsatsanis and G.B Giannakis, "Object and texture classification using higher-order statistics", IEEE Trans. Acoust Speech Signal Process., Vol. 38. July 1990, pp. 1284-1296.