Bispectral techniques for spherical functions

zyx zyxwvuts zyxwvutsr BISPECTRAL TECHNIQUES FOR SPHERICAL FUNCTIONS Ra,makrzshna Kakamla', Bruce M. Bennett", Geoflrey J . Iverson", and Mzchael D 'Zmura** 'Dept. of Elect. & Electronic Eng., Univ. of Auckland, Private Bag, Auckland, New Zealand **Depts. of Math. and Cognitive Sci., Univ. of California, Irvine CA 92717, USA zyxwv zyxw ABSTRACT 2. DEFINITIONS Spherical functions arise in geophysics, medical imaging, and computer graphics. This paper addresses two problems involving spherical functions: determining when two spherical functions are 3-D rotated copies of each other. and averaging several noisy observations of a rotating spherical function. Our solution t,o both problems uses the spherical bispectrum, which is the generalization of the wellknown Euclidean bispectrum. In this paper, we formulate the spherical bispectrum and show t h a t it is invariant under 3-D rotation of the underlying function and unbiased in the presense of additive Gaussian noise. We demonstrate an algorithm for recovering spherical functions from their bispectra. T h e bispectrum of any function f on R" is the Fourier transform of t h a t function's triple correlation. In this section, we demonstrate how to obtain the spherical triple correlation from its well-known analogue on Euclidean space. zyxwvut zyx zyxwvutsrqpon A . Euclidean triple correlation Let f be any complex-valued function in Ll(R"). Its triple correlation a3.j is defined by the integral a 3 , f ( s ,1 ) = J,. f ( z ) * f ( z + s)f(. + t)dz, (1) where the asterisk * denotes complex-conjugation. Clearly, the triple correlation is invariant under translation of f , i.e., if there exists t such that h ( z ) = f ( z t ) for all z, then t i 3 , / , = a 3 , f . T h e bispectrvm - 4 3 , f is the Fourier transform of the triple correlation. One finds easily t h a t + 1. INTRODUCTION Spherical functions (functions whose domain is a sphere) are useful for modelling the earth's gravitat,ional field [I], describing the deformation of the heart surface during a cardiac cycle [2], and modelling reflectance from a spatially extended light source [3]. There are two general problems involving spherical functions t h a t we address in this paper: how to determine when two given spherical functions are simply three-dimensional (3-D) rotated copies of each other; and how to average several noisy observations of a single spherical function t h a t is undergoing an unknown 3D rotation in between successive observations. For similar problems, when the functions involved are defined on the real line, bispectral techniques are appropriate for the following reasons [4]: (1) the bispectrum of any deterministic function on the real line is invariant under translation of the function; (3) the expected bispectrum of any zero-mean Gaussian noise process is identically zero; (3) the bispectrum of most deterministic functions on the real line contains enough information to enable recovery of the underlying function. In this paper we formulate the appropriate bispectrum for spherical functions, show that the new spherical bispectrum possesses attractive properties similar to (1)-(3) above, and demonstrate bispectral applications to the matching and averaging tasks mentioned above. '43,f(U, U) + U)*, = F(u)F(v)F(ffl (2) where F is the Fourier transform of f [4]. In many cases, the bispectrum A3,f is unique to the underlying function f and its translates, i.e., if for some function h we have that A3,h = A 3 , j , then there exists 1 such that h ( z ) = f ( z 1) for all 1. For example, this is true if f is nonnegative and has compact support; for a detailed review of other cases where uniqueness holds. see [SI. The triple correlation for deterministic functions, as calculated in e q . ( l ) , is unbiased in additive Gaussian noise. To make that precise, suppose t h a t we observe over a compact region K of IR" a signal r of the form r(z) = f(r) n ( z ) , where f is a deterministic function and n is a sample from a zero-mean stationary Gaussian process. If we subtract from r its mean value 7 over Ii, and set r' = r - 5,then we find that E [ a s , + ( s :t ) ] = a3,fl(s11). + zyxwvutsrqpo zyxwvu This work was supported by NSF Grant DIR901.1278 to the Jnstitute for Mathematical Behavioral Sciences, University of California, hvine, R. D. Luce, Director. *E-mail address: R.KAI(ARALA~CCUl..~UI<UNI..~C.~Z + Here f ' denotes the zero-mean function f - 7. B. Spherical triple correlation T h e Euclidean triple correlation in eq.(l) is the integral of the Euclidean function f multiplied by two independently translated copies of itself. Therefore, the appropriate spherical triple correlation should integrate the spherical function f multiplied by two independently rotated copies of itself. To make _the analogy even more direct, we define a new function f on the 3-D rotation group S 0 ( 3 ) , in such a way that if f undergoes a 3-D rotation, then f undergoes a corresponding left translation on SO(3). We show later how zyxwvuts IV-216 0-7803-0946-4193$3.00 0 1993 IEEE - ~~ zyxwvutsrq zyxwvutsrqpon zyxwvutsrqp zyxwvuts zyxwvutsrq zyxwvut zyxwvutsrqpo zyxwvut zyxwvut to const,ruct f from f, but for now let. f be an arbitrary function on SO(3). We define the triple correlation o f f to be the integral G,,j(S,T) = J ~ ( R ) * ~ ( R s ) ~ ( R T ) (~3R) . SO(3) Here dR is the unique normalized Haar integral on S 0 ( 3 ) , which is invariant under either left or right translation of f, i.e., J f ( R R o ) d =~ J J ( R ) ~ =R J j ( R ~~ R ~. Here, the dagger t denotes matrix adjoint. T h e translation property of th_e Fourier coefficients is important in the bequel: h ( R ) = f ( S R ) for all R if and only if H ( ! ) = F ( t ) D e ( S )for all The triple correlation is a function on SO(3) x S 0 ( 3 ) , and therefore requires for its Fourier expansion Kronecker products of the Dt matrices. T h e (k,C)-th Fourier coefficient matrix of u 3 , j , which we denote A , , j ( k , C ) , is computed by the integral e. J J (4) (The explicit form of the Haar integral in terms of Euler angles is contained in standard text,s [6].) By applying (4) to eq.(3), i t is easily seen that u 3 , i is invariant under left translation of f on SO(^). We now apply eq.(3) to spherical functions. Let S2 denote the sphere of unit. radius in IR3, and let h = [U. 0, l]‘ denote the North pole (1 denot.es transpose). To each function f on S’, let denote the function on SO(3) defined by the rule j ( R ) = f ( R h ) . We define the spherical triple correlation of f to be u 3 , j in eq. (3) for the corresponding function f on SO(3). This triple correlation is invariant under 3-D rotation of f. To see that, suppose h is another spherical function that is a 3-D rotated copy of f. Then there exists a rotation Ro such that h ( z ) = f ( R o z ) for all z E S’. Thus k ( R ) = h ( R k ) = f ( R 0 R A ) = f ( R o R ) ,or equivalently, the two functions are left trauslates of each other. and consequently u3,,j = u 3 , , f . Now suppose we observe r ( x ) = f ( x ) + n ( x ) on S 2 , where f is a determinstic function and .n is a sample function from a Gaussian, zero-mean, stationary process N . As above, i t may be shown that a , , j ( S , T ) [Dk(S)t @ De(T)’] dSdT, SO(3) SO(3) where @ denotes the matrix Kronecker product. By the Clebsh-Gordan formula, the matrix D k ( S ) @ D e ( T reduces, ) when S = T , into the following direct sum of smaller matrices : C k e [Dk+e(s)@ Dk+e-i(S) @ . . . @ Dp-cl(S)] cLe ( 5 ) Here C ~isCthe unitary matrix of Clebsh-Gordan coefficients (independent of S ) , and & denotes matrix direct sum. The discussion above gives us the machinery required to compute the spherical bispectral coefficients [7, pg 781. Proposition 1 Let f be any Lz functzon on SO(3) wuath Fourier coeficzents E ( Q , C 2 0. For each k und e, we have c) .J?,J~. = [ E ( k )a E ( [ ) ]ckC [ E ( k + e)t@ P ( k + e - I)+ . . . e E(p -el)+] zyxwvutsrq E [~3,i.’(S,T ) ] ~ 3 , j , ( ST), , with i‘ and f’ denoting respectively the extensions to SO(3) of the mean-subtracted functions 7 - F and f -7, the mean being evaluated on the sphere [7, pp 21-22]. 3. BISPECTRUM FORMULA In this section, we calculate the Fourier represeiitation of the spherical triple correlation. Fourier analysis on the group SO(3) is determined by the irreducible unitary r e p resentations of the group [ 6 , Ch. I]. T h e latter are matrix valued functions De. e 2 0 , whose coefficients are called generalized spherical harmonics. For each 1, the function Dl maps SO(3) into the group of unitary matrices of dimension ?C 1, in such a way that De(RS) = De(R)De(S) for all S, R. For any Lz function f on SO(3), we have the Fourier series expansion + CO e=o / 4. BISPECTRUM RECOVERY ALGORITHM The formula in Proposition 1 suggests that i t is possible to extract the underlying Fourier coefficients recursively from the bispectrum. We demonstfate below an algorithm accomplishing this for functions f defined on S 0 ( 3 ) , and later. we show how the technique applies to spherical functions. Our algorithm makes use of the following facts from matrix theory ([a]). First. any positive definite matrix H has unique “positive square root”, i.e., a positive definite zyxwvutsrqpon so(3 ) Hi 1 1 1 where tr denotes matrix trace, and the Fourier coefficient matrices F are obtained from the matrix-valued integral: F(e)= Intuitively, the spherical bispectrum formula follows from the Euclidean formula (2) by taking into accout the Kroneckerproduct decomposition formula (5) on the group SO(3). (For the group R, the Kronecker-product decomposition of irreducible unitary representations is simply the familiar expression e*xz.etu2= 1 If f is any function on the sphere, and f is its extension to SO(3) by the North pole method of $11, then we call A 3 , y the .spherical bispectrum of f. matrix such that H f H f = H . Second, any nonsingular matrix A has a unique polar decomposition A = H+I:, t r [ f ( C ) D e ( R ),] f ( R )= cle. f(R)De(R)‘dR. where H+ = (AA’):. and U is a unitary matrix. Now let L > 0, and let f be any real-valued_ function on SO(3) whose Fourier coefficients are such t h a t P ( I ) is a nonsingular matrix for each e L . and furthermore, P ( e ) = 0 if e > L. Since f is real-valued, it is easily shown that < IV-217 zyxwvutsrqponm zyxwvutsrqponm zyxwvutsrqpon zyxwvutsrqpo zyxwvut zyxwvutsrqpo zyxwvutsr det [P(1)] is a real number. Assume, for now, that it is positive. Our algorithm for recovering from its bispectrum A 3 , j proceeds in three steps. 1. Since f is real-valued, F ( 0 ) is a real number, and by Prop. 1, A 3 , j ( O , 0 ) = P ( 0 ) 3 . and thus E ( 0 ) = y.- 2 . We estimate the first Fourier coefficient matrix E( 1). Because E ( 0 ) is. by assumption, a nonzero real number, we use the formula in Prop. 1 and find t h a t as "side information" along with the triple correlation, then we obtain a complete (left) translation-invariant description for any real-valued bandlimited function on S O ( 3 ) . Note that det[E( l)] remains invariant under translation on S 0 ( 3 ) , i.e., if f ( R ) = i ( S R ) , then B ( 1 ) = h ( l ) D l ( S ) , but since d e t [ D ~ ( S ) ]= $1 ([9, pg 471) we obtain that det[P(1)] = det[b(l)]. T o sum up, any real-valued bandlimited function f on S 0 ( 3 ) , whose coefficient matrices are all nonsingular up to the bandlimit, can be recovered completely-up to a single left translation on S0(3)-if both its triple correlation and the value of det[E(l)] is known, and the algorithm described above is used. 5. SPHERICAL FUNCTIONS T h e matrix on the right hand side above is posi.4, f ( l 0 ) 1 tive definite. Let F ( l ) = be its positive square-root. It can be shown ([7, p p 116-1181) that if E(1) is constructed in this way, then there exists S in SO(3) such that (w): P(1) = P ( l ) D I ( S ) . 3. If L = 1, then we are done. Otherwise. t!ie following recursion produces matrices F ( 2 ) .. . . , F ( L ) , such that for all 2 5 t' 5 L we have t h a t p(1)= E ( c ) D e ( S ) for t h e same S as in Step 2 . Since we know p(l) ant1 i13.j(1, I ) , we obtain p ( 2 ) from the upper-left 5 x 5 submatrix of the following 9 x 9 matrix: Explanation: all terms above are known, and if we substitute for F(1) and A 3 , j ( 1 ,l ) , and use the reduction formula (5), then we find that The algorithm described in the previous section is formulated for functions on S 0 ( 3 ) , but we need additional steps to use it for spherical functions. Let f be the extension t o SO(3) of a spherical function f via the North pole mapping. The Fourier coefficient matrices E ( O ) , E ( 1 ) , . . . of any function f obtained in chis way have at most rank 1. T o prove this, consider t h a t f is invariant under any rotation Q t h a t leaves the North pole fixed, i.e., ~ ( R Q=) ~ ( R Q & =)~ (RL= ) J(R). The set of all rotations Q t h a t fix the North pole forms a subgroup Q of SO(3). T h u s we can average each Fourier coefficient matrix E ( [ )over Q in without changing the value of the coefficient. It is easily shown t h a t this implies E ( [ ) = PtE(t'),where Pe = J, De(Q)dQ. zyxwvutsrqpo [ D z ( S ) ' E ( 2 ) ' ]@ [ D , ( S ) t E ( l ) ' ] E(0) T h e upper left 5 x 5 submatrix of the right hand side is exactly the matrix [ F ( ? ) & ( S ) ] and we set its adjoint equal to p ( 2 ) . Having obtained P ( 2 ) , we estimate $'(e) for any e > '7 from the adjoint of the upper left ( 2 l + 1 ) x ('7e+l) submatrix of the following matrix: ', T h e same argument as above shows t h a t p ( l ) = P ( t ) D t ( S ) .After iterating the recursion until e = L , the function f on SO(3) obtained by Fourier series expansion with the coefficients E(()),p ( I ) , . . ., p ( L ) , is such t h a t f ( R ) = f ( S R ) for all R. T h e assumption t h a t det[F(l)] > 0 is not critical. Instead of selecting P ( 1 ) t,o be the positive definite square root of F ( l ) P ( l ) ' , we may choose P ( 1 ) to be any square root such that d e t [ P ( I ) ] = det[F( l)],e.g., by multiplying the top row of the positive definite square root matrix by -1 if necessary. We do not know det[F(1)] a priori, but if we store it The matrix Pe is an orthogonal projection, whose elements are zero everywhere except for the exact center-of the matrix, which has the value l ; thus t h e rank of F(E) cannot exceed one [7, pg 951. Although _the Fourier coefficients of an extended spherical function f are singular, we may augment them so t h a t they become nonsingular, compute the resulting bispectrum, and then use the recovery algorithm in fIV. This allows us to average multiple observations of a rotating spherical function, a procedure t h a t is described in $VI. We focus now on the method of constructing nonsingular Fourier matrices. Since F ( t ) = PeE(t), the entries of the 2t+l-dimensional matrix E ( t ) are zero except for the middle row. Assume that at least one of the middle row elements is nonzero; then we may substitute for t h e 2 l remaining rows linearly independent vectors so t h a t the resulting matrix, denoted P(t'),is nonsingular. (We are finding a nonsingular matrix P [ e ) such t h a t $ ( l )= P ~ p ( l ) . Moreover, ) -we may choose the additional vectors so t h a t t h e function f defined by the coefficients P(O), P ( l ) , - . . ., P ( L ) is real-valued. Now the function f is uniquely determined by its bispectrum A 3 , j up to a 3-D rotation, i..e, from A3,j we can recover Fourier coefficients P such that for spme S a n d all 0 5 l < L , we have P ( l ) = p ( l ) D t ( S ) .We then obtain the IV-218 zyxwvutsrqpon zyxwvutsrqpon zyxwvutsr zyxwvutsrq zyxwvuts zyxwvutsrqp zyxwvutsrqp zyxwvutsrqp original Fourier coefficients (up to a 3-D rotation) by composition : P e p ( [ ) = f i ( t ) D t ( S ) . Thus the bispectrum of f uniquely determines and consequently also the original spherical function f , up to a 3-D rotation. i. 6. APPLICATIONS We outline two applications of the spherical bispectrum, one to matching spherical functions, and the other to averaging multiple views of a rotating Spherical function. Since the spherical bispectrum is invariant under 3-D rotation of the underlying function, we check whether two spherical functions y and q are simply 3-D rotations of each other by comparing their bispectra. T h e bispectrnm test has two appealing properties. First, the recovery algorithm in $4, although formulated for functions on SO(3) with nonsingular coefficients, shows t h a t bispectral coefficients are potentially powerful enough to uniquely determine the underlying function and its translates. Second, the insensitivity of bispectral coefficients to additive Gaussian noise suggests that any test involving matching bispectral coefficients should be insensitive to measurement noise. In practice, however, the bispectra of rotated copies of the same fuiiction will not match exactly, and thus we need to check if their difference is less than some threshold, where the latter is determined by measuring or calculating the probability distribution of the bispectral coefficients. We leave this topic for future investigation. Suppose now t h a t we have several noisy observations fi . . ., f n , of a single spherical function f. where the function may be rotating by an unknown amount in between successive observations, and the noise is independent in each observation. Since the signal is rotating, simply averaging the observations would average out the signal. Instead, we may use the following bispectral technique to average the noise without averaging out the signal. For each observation, we compute the spherical bispectrum, and average the resulting bispectra over all the observations. In doing this, we are averaging noise without averaging out information about the signal. Unfortunately, we cannot apply the algorithm of $4 to recover the original signal from the averaged spherical bispectrum < A 3 , i >, because the underlying Fourier coefficient matrices are singular. However, we may recover the underlying signal if we select any one of the observations, say fn, form the augmented function f n on SO(3) with nonsingular coefficients, and compute the bispectrum A,,fn. T h e nonzero coefficients of < A 3 , j > form a proper subset of the coefficients of A,,fn, and thus we may improve signal-to-noise ratio in the reconstruction by substituting in the appropriate locations of each matrix A,,,* (k,e) the coefficients of the averaged matrix < A , , j ( k , t) >. Numerical results obtained by using this procedure will be reported in future work. ing functions on the 3-D rotation group from their bispectra, and we describe how that algorithm may be used with spherical functions by augmenting the Fourier coefficient matrices. Finally. we describe applications of the spherical bispectrum to the problem of checking whether two spherical functions are simply 3-D rotations of each other, and to the problem of averaging multiple views of a single rotating spherical function. 8 . REFERENCES [l] R. A. Nash and S. I<. Jordan, “Statistical geodesy: an engineering perspective,” Proceedings of the IEEE, Vol. 66, No. 5, pp 532-550, May 1978. [2] C. W. Chen and T. S. Huang, “Left ventricle motion analysis by hierarchical decomposition,” in Proceedings ICASSP-92. Vol. 3. pp 373-?76, 1992. [3] &I. D’Zmura, “Shading ambiguity: Reflectance and illumination,” in Compututional models of visual processing, M. S. Land? and A. J. Movshon, Eds., Cambridge: MIT Press, pp 187-207,1991. [4] A. W. Lohmann and B. Wirnitzer, “Triple correla- tions,” Proceedings of the IEEE, Vol. 72, p p 889-901, 1984. zyxwvutsr . [5] J. I. Yellott Jr., and G. J. Iverson.“Uniqueness results for generalized autocorrelations,” Journal of the Optical Society of .America -4,vol. 9, p p 388-404, 1993. M. Gel’fand, R. A. Xinlos, and Z. Ya. Shapiro, Representations of the rotation and Lorentz groups and their applications, Oxford: Pergammon Press, 1963. [6] I. [7] R. Kakarala, Triple correlation on groups, Ph.D. The- sis, Department of Mathematics, University of California, Irvine, 1992. [8] P. Lancaster and M. Tismenetsky, The theory of matrices, Orland: Academic Press, 1985. [9] L. C. Biedenharn and J. D. Louck, Angular momentum in quantum physics: Theory and application, Reading, MA: Addison-Wesley, 1981. zyxwvutsrqp 7 . SUMMARY In this paper, we formulate the appropriate bispectrum for spherical functions, and show that it in unbiased in Gaussian noise and invariant under 3-D rotation of the underlying function. We demonstrate an algorithm for recover- IV-219