Digital restoration of images in the presence of correlated noise

Comput. & Elect. Engng,

Vol. 3, pp. 65--74. Pergamon Press, 1976. Printed in Great Britain

DIGITAL RESTORATION OF IMAGES IN THE PRESENCE OF CORRELATED NOISE B. P. AGRAWAL General Dynamics(ElectronicsDivision),Orlando, FL 32802,U.S.A.

and V. K. JAIN Departmentof Electrical and Electronic Systems,Universityof South Florida, Tampa,FL 33620,U.S.A.

(Received 8 May 1975) Abstract--Digital restoration of images, basically a problem in two-dimensional deconvolution, is complicated by the large dimensionality and inherent ill-conditioning of the describing equations. In this paper we develop a constrained-least-squares solution assuming that the degrading phenomenon can be modeled as a linear isoplanar system, and the additive noise as a Pth order Markov process. To achieve a practical computer algorithm, the block-Toeplitz and block-circulant properties of the imaging equation matrices are exploited. Specifically,these properties are utilized to adapt Kutikov and Akaike's algorithms for the solutionof the equations.Not onlyis the modelmoregeneralthanhas heretoforebeenconsidered,but our approach results in considerable savings in memory storage compared to the existing approaches. 1. INTRODUCTION Potentially, one of the most rewarding areas of image processing by digital computers is the restoration of degraded images. The objective therein is to undo the loss in quality suffered due to a variety of causes, such as defocusing, atmospheric turbulence, and certain linear motions of the source image during photographic exposure. Other possible factors are: motion of the optical medium, e.g. sea water, and backscattering media such as smoke and haze. Several deconvolution methods have therefore been proposed[I--6] for inverting the degradation. A linear isoplanar degradation is generally assumed wherein the additive noise is considered as an uncorrelated process, or is ignored. Due to ill-conditioning of the deconvolution equations, even a small amount of noise creates high-frequency oscillations of significant amplitude in the solution image. Hunt [1] has therefore attempted a constrained-least-squares solution and has provided a detailed, successful application procedure. The present work provides a solution to the image restoration problem in a more general setting where we permit the additive noise to be a Pth order Markov field [7]. Furthermore, our algorithm exploits the block-Toeplitz and block-circulant properties of the imaging equation matrices so that memory space requirements are reduced by a factor of 10 or more for the computation of the restored image. This is made possible by adapting Kutikov's and Akaike's methods for finding the inverse of matrices with such special properties.

Digital image restoration problem For isoplanatic degradations such as those mentioned above, or those introduced deliberately for bandwidth compression during transmission[17], a linear imaging model is

g(k,l)=

~

h(k-m,l-n)f(m,n)+v(k,l).

(1)

m,n = -~

Here, [(k, 1) is the source image, h(k, I) is the point spread of the image formation system, and g(k, l) is the received or degraded image. We assume, without loss of generality, that the source image and the point spread of the image formation system are defined over the square lattices L.a and Ib,b respectively and are identically zero outside. The additive noise v(k, l) is a Pth order Gaussian Markov random field[7] with

Ev(k, I) = 0 Ev(k, l)v(v, I~) = I q([k - v[, [1 - / z l) for ( k - v)2+(I-/~)2
(2) (3)

66

B . P . AGRAWAL and V. K. J~lN

The digital image restoration problem posed here is: Given the degraded image g, the point spread matrix h, and the correlation function q, estimate the source image f from the model (I). 2. V E C T O R - M A T R I X

FORMULATION

In order that high-speed convolution algorithms may be applied, the matrices f, h and k, are augmented with zeros so that each is defined on an extended lattice l ~ . N is chosen equal to or greater than a + b - I to avoid wrap around during convolution. Stated explicity, let the matrices f, h and g be redefined as t

°I1 f(k,I)

Z

0

0 h=h

h (k I)

o

0

0

I

I N-I

I E

J

I

+

I

The augmented noise matrix v is v =[v(k,l)], k,I =0,1 . . . . . N - I . The discrete twodimensional convolution in (1) can then be replaced by circular convolution and the model becomes N-I

g(k,l)= ~, h ( ( k - m , l - n ) ) f ( m , n ) + v ( k , t ) . mm

(4)

-0

Equation (4) can be rewritten in a more desirable, vector-matrix form as follows: Step 1. By ordering the elements of matrix f lexicographically, construct an extended column vector F of length N 2. That is, let F' = [fi,, f't . . . . . /L ,] where/k is the (k + 1)th column vector of the augmented matrix f(k, I). The primes denote transpose operation. Step 2. Similarly, construct the extended vectors G and V from matrices g and v. Step 3. Form a matrix H, of size N 2 x N 2, with N 2 partitions or blocks. Each block contains N x N entries• The blocks are ordered as follows:

YG H =

.

1

~N

(51

2

Each partition of H is generated from a single column of the matrix t7. Specifically, h(O,j)

I

h(1,j).

=

Lh(N

h(N-l,j) h(O,j)

. .h(1,j)

• .h(2,j) (6)

- I, y)

h(N

-,~ j). . .h(O,j)

Observe that ~ is a circulant matrix• That is, each of its rows is a right shift of the preceding row, except the first which is a right shift of the last row. Also, H is a circulant matrix in terms of its blocks•

Digital restorationof imagesin the presenceof correlated noise

67

The preceding construction of the vectors F, G and V, and the matrix H enables us to rewrite (4) in a vector-matrix equation of the form G = ~'F + V.

(7)

The noise vector is characterized by the statistical properties

EV= 0 EVV'

=

(8)

Q.

(9)

The correlation matrix Q is a highly structured matrix. In fact, it has the form

Q, Q=

Qo...QN-2 .

. . . .

(lo)

. . .

Qo ~ QiN, QN-2 Q' ... QN-,[ where each block Qk, k = 0, 1. . . . . N - 1 is an N × N matrix. Their entries represent the correlation between the components of the noise vector V and can be determined by use of (3). 3. PROPERTIES OF IMAGING EQUATION MATRICES The correlation matrix Q and the transmission matrix H carry the bulk of information about the imaging system represented by (7). They both possess some very interesting properties which will help facilitate our computation of the restored image.

Correlationmatrix It is clear from (10) that Q is an Nth order block-Toeplitz matrixL In addition it exhibits the following properties: (i) Each block of the correlation matrix Q is itself a Toeplitz matrix. (ii) Q is symmetric both element-wise and block-wise. Further, Q is a block-persymmetric matrix$. As a direct consequence of these properties Q-~ is also symmetric both element-wise and block-wise, and is block-persymmetric. (iii) An immediate corollary of (ii) is that both Q and Q-' are completely described by their respective first rows (or first columns). (iv) Finally, being a correlation matrix, Q is a positive definite matrix as is also Q-'.

Image TransmissionMatrix H, and the product H'H The matrix H displays the following properties: 6) H is an N × N block-circulant§ matrix of block-size N x N, i.e. each block-row of matrix H is obtained by a circular right shift of the block-row just above it and the first block-row is the circular right shift of the last block-row. (ii) Each block of the matrix H is itself a circulant matrix. (iii) Property (i) implies that H is a block-Toeplitz matrix. Properties (i) and (ii) together imply that H is completely specified by its first row (or column). (iv) Since the product of two block-circulant matrices is also a block-circulant matrix, the product H'H is a block-circulant matrix. Further, note that each block of H'H is a circulant matrix and, therefore, H'H is completely specified by its first row (or column). (v) Finally, it is observed that the inverse matrices H ' and (H'H) ' are block-circulant

~An N × N block-matrixP.whose i,./block,P~,is a functionof (i -/') only,is calledan Nth orderblock-Toeplitzmatrix. :~Amatrix symmetricabout its cross diagonal is said to be persymmetric. Block-persymmetry,likewise, means that the matrix is block-symmetricabout the cross block-diagonal. §Each block row of a block-circulantmatrix can be obtained from the precedingblock row by one block shift to the right; one exception is the first block row which can be obtained from the last by one block shift to the right. Note that a block-circulant matrix is always block-Toeplitz,although the converse is not necessarily true.

68

B. P, AGRAWAL and V. K. JAIN

matrices and each of their blocks is also circulant. Therefore, both H ' and (H'H) ' can be described by their respective first rows (or first columns).

4. CONSTRAINED-LEAST-SQUARESESTIMATION We now return to the main problem of interest, image restoration. An obvious method of estimating the source image is to solve (7) by least-squares method. Like all direct deconvolutions, however, this solution is unreliable. Because of the inherent ill-conditioning of the matrix H[I 1, 12] such solutions are very sensitive to perturbations in the image G, and are known to result in undesirable high-frequency oscillations in the estimated image. Hunt[l], following the idea suggested by Phillips [8], incorporated a smoothness measure for the estimated image. He also reduced the time of computations by an ingenious use of the fast Fourier algorithm (FFT) to carry out certain matrix inversions. The success achieved in this work, therefore, represents a significant advance in the practical solution of the image restoration problem. The present development provides a different solution. The advantages realized are: (a) by permitting the noise in the model to be a Pth order Markov field the formulation is more realistic (Hunt considered uncorrelated noise in his formulation: his solution is inapplicable to the present general case), and (b) our memory storage requirements are of the order of N ~ for an N × N image resolution (rather than N ' as in Hunt's solution). This latter fact makes it possible to carry out image restoration on moderately large computers. Constrained-least-squares formulation From eqn (7) it follows that (G-HF)'Q '(G-HF)=

V'Q i V

Ill)

where we note EV'Q 'V = N:.

(12)

It appears reasonable, therefore, to require that the estimate f" satisfy the constraint

( G - H F ) ' Q ' ( G - H F ) : N ~.

113)

To quantify the smoothness, actually the lack of it, we adopt a measure originally proposed by Phillips, namely the inner product {Ct+, CF}.

t14i

Here, C is an N 2 × N~" roughness matrix with an approximate second difference structure. This matrix can be constructed from the N × N matrix

0 (~l

0 O

0

0

~,

(15/

O

where C1 is the discrete Laplacian operator. Now, convolution of an N × N image matrix, say ], with CLe is equivalent to matrix multiplication of the corresponding lexicographic vector, F, with

~ '., C=

f

O

i' (:" .~ Ci

... (il2t

hO o :[: c',;

•

(16)

Digitalrestorationof imagesin the presenceof correlatednoise

69

The N × N block C. is generated by circulating the first column of CLe,i.e. it has a~ as the entries along the lower sub-diagonal, an a, in the upper right hand corner, and zeros elsewhere. Similarly, C, is generated by circulating the second column of CL,, and C2 = Co. The numbers a, and a2 would have to be 1 and 4 respectively in order that the 3 × 3 matrix CL be the discrete Laplacian operator. In practice, however, this choice must be avoided. Slightly perturbed values must be used, for otherwise a singularity is produced in the solution. Hence the remark above that C has an approximate second difference structure. Problem statement. Given G, H, and Q, compute the estimate/6 which minimizes the inner product (14) subject to the constraint (13). Solution. A straightforward Lagrangian minimization yields /~=(H'Q ' H + A C ' C ) ' H ' Q 'G.

(17)

In principle the problem is solved by eliminating A from (17) with the help of (13). However, the matrix inversion in (17) presents a serious computational difficulty which we will attempt to resolve in the subsequent sections. First, we modify (17) by use of the matrix identity

(A +B'DB)-' = A - ' - A - ' B ' ( B A 'B'+D ~) 'BA-',

(18)

where A and D are arbitrary nonsingular square matrices. The use of this identity in (17) yields i~"= [S-SH'(HSH' +Q) 'HS]H'Q-' G

(19)

S ~= AC'C.

(20)

where

At first glance the right hand side of (19) seems to demand formidable memory storage capacity and computation time. However, the matrices Q, S and (HSH'+Q) are all symmetric block-Toeplitz matrices and, therefore, their inversion reduces to computing the first block-row of their respective inverse matrices. As for memory storage, we recognize that S, HSH' and Q are completely specified by their respective first rows (or columns), hence only these need be stored. To compute the inverse of the matrices Q, S and (HSH' + Q), we shall adapt the block-Toeplitz matrix inversion algorithm due to Kutikov[13] and Akaike[14].

5. INVERSION ALGORITHM FOR SYMMETRIC BLOCK-TOEPLITZ MATRICES For the discussion below, the following conventions shall be observed. Lower case symbols Denote a block of size d x d. Block size will be kept fixed in the entire discussion. Upper case symbols Denote a block-column. Supernumeral in parentheses will indicate the number of blocks in the column. Denote a square block-matrix. A supernumeral in parentheses will Bold face upper indicate the order of the block-matrix. case symbols Denotes block-transposition, i.e. replacement of ij block by fi block in a block-matrix. Denotes the transposition of elements of a matrix. Denotes the reversal of ordering (of blocks in a block-row or block-column.) A symmetric (elementwise) (N + 1)st order block-Toeplitz matrix L'N+'> and its inverse A'N+'~ can be represented as follows:

~lo i L'(N)~

(21)

B. P. AGRAWM. and V K JatN

70

and ¢{N)

[L ,N~,,] , = A,N~ ,,= [..a.~2.L.i.O...... ] I A'N'!B'~'I .

.

.

.

('}~}

Note that 1o is a d x d matrix. 1-'N' = [1~.1~ . . . . . IN] {the block transpose of 1 .....}. is a I :y N block-matrix and L 'N' is an Nth order block-matrix. Analogously, the block-transpose of L ~N~'' and its inverse U '~ .... can be represented as:

"

~{~r, ]

{23}

: [ u,N)7 v,N, 1.

~24}

=[I~,,NJj

and

Note that the operations of block transposition and matrix inversion do not commute (unlike the elementwise transposition and matrix inversion, which commute). Recursive c o mpu t a t i o n of the first block column of A ...... As mentioned earlier Akaike presented an iterative algorithm[14] for the computation of the inverse of a block-Toeplitz matrix. To summarize it. let •'~2' ..... = l u . , , , ]

'.4 ..... '~

{25}

UY ' + ' , = [u,,+,] 'U .......

{26)

and

Using this notation, it can be shown that aN +,. A ,N,, UN + ~ a n d U ~ ' can be generated recursively as follows: t{n

UI,'N~"=IU,,

P

t(n

;0I-(U,.

) --{rt }

L

-- t { n

+l,,~,}a,,[A.

'iil

{27)

ft~, (°*', =[01A. ' - " " ' ] - ( A -. ' { " ' -L" " '+ I,, , ,}u,,IilU. ' .... I

128)

' = Iu,, I ' -{ U,', (..... } , , . , ( A 2 ' " ' " ) , l u , , 1

~29)

lu,,,,]

and [a,,,,] ' = [ a , , ] ' - ( A . --

p{n ~ I)

),{U v.{.n ,. . I } )

z

,(a.

l

(30)

where the suffixes 1 and n +1 denote the first and the last blocks of U,,''~''' and A,,''''~' respectively. An identity block is denoted as i. The recursive relations (27}-(30} are initialized by considering n = 1 and n + I = 2:

1/'' l;1 l/,

{Bt

l,,l

and A,~,

[

(l,,-t;to

[-(t,,-t,l.

't,) '

- 1,, 'l'I(lo-

1,1,, '1;} '1.

(1o- l,li;'l'D

'8)'l,l,,'

(32}

Hence, [a21' = (l,, - 1,1{, ' '1,) A{"=(I,,-

I,I{, 'l't) '1,1,,'

[u21 ' = (l,~-l,l,, '1;t '

(33) {34t {35

Digital restoration of images in the presence of correlated noise

71

and (36)

U ' ' ' = (lo- l'mlo-ll,)-'l'do-I.

Using these initial values and eqns (27)-(30) one can recursively compute the first block-column of A~N+'.

Computation o[ remaining blocks of A~N÷" Kutikov and Akaike have shown that A t kN++l , j,+ l

_ AtN~ ±I ~~ zal t N ~1a1 a'tN) t'Xkd

--

--

UtN)U,u~N~)k,j

k, j = l, 2 . . . . . N.

(37)

With the help of (37) and the fact that A"~÷" is element-wise symmetric and block-persymmetric, all of the subblock-matrices can be computed easily. 6. A PRACTICAL SCHEME FOR CONSTRAINED-LEAST-SQUARES ESTIMATION The constrained-least-squares solution to the image restoration problem was shown to consist in solving simultaneously, /~ = [S - SH'(HSH' + Q)-IHS]H'Q-t G

(19)

and the constraint equation in (13), rewritten as ~b(,~) = (G - H~e)'Q-'(G - HF) - N 2 = 0.

(13')

An explicit solution of the above pair does not seem possible; therefore we use Newton-Raphson iteration to find A: A.,+,

= ,~., - ,/,(x,,)/,//(,A,,).

(38)

Here, Ak represents the kth step approximation to the true root of $(A). The iteration can be started with Ao> 0 for k = 0 and continued until a desired degree of accuracy is reached. A unique solution exists if G'QG > N 2. This can be demonstrated easily by noting that ~b(0) = - N 2, and by showing that (d~b/dA)> 0[15]. The details are left to the reader.

Steps in iterative algorithm Step 1. From the transmission matrix H, the second difference matrix C, and the correlation matrix Q. Compute the inverse of Q using the algorithm for the inversion of symmetric block-Toeplitz matrices. Step 2. To start the iterative algorithm, choose a positive initial value for A, say Ao. Step 3. Compute (B'SH+Q) -t using the algorithm for the inversion of symmetric block-Toeplitz matrices described in Section 5. Step 4. Compute the constrained-least-squares estimate in accordance with (19). Step 5. Compute ~b(A) and ~'(A). If (a) [~(A)[ > a, increment (or decrement) A according to (38). The constant a is preset and is a measure of accuracy with which the constraint (13') is to be satisfied. Return to Step 3. (b) [~(A) I < a, stop the iteration. 7. COMPUTER SIMULATION A FORTRAN IV program was written for the IBM 360/65 computer to test the iterative algorithm proposed in the previous sections. First, several tutorial examples with low image dimensionality were run; one such example is given in [15] together with all the intermediate results. Here we present the results for a 24 x 24 image which required 90 K words (360 K bytes) of core memory. The first iteration of the program, steps 1-5 of the algorithm, was completed in 4450 sec. The successive iterations took approx. 900 sec each.

72

B . P . AGRAWALand V, K. JAIN

The object image and the point spread of the imaging system chosen ;ire:

f(k,/)

10 for 20 for =

0-<-k-<10 and 10
20 for 0-
and

()~/fl0 10
10<1¢--.20

and

h(k 1)= {0.0~

for O<-k,l<-5 elsewhere.

The additive noise field v(k, I) is modeled by a first-order homogeneous Gauss-Markov random field with

Ev(k,I)=O and l

t

for k = v . l = #

Et,(k,I)u(u,~)= 0-1 for O<(k-v)2+(l-bt)?<.l I()

otherwise.

The actual numbers having the above characteristics were generated by Woods' method [7]. The degraded image matrix g(k, l) was obtained by adding the noise matrix to the circular convolution of the object image with the point spread matrix. The source image and the degraded image are shown in Figs. I and 2 respectively. The accuracy sensitivity parameter o~, used in Step 5 of the algorithm, bad been set to reflect 5% error in the satisfaction of the constraint, i.e. oL= 28.8. For A = 0.0102 the value of O(A) was computed to be 563.8. The restored image matrix / as computed by the program is shown in Fig. 3. For numerical comparison the degraded image matrix and the restored image matrix are given in Table 1. 8. CONCI+USIONS A computer algorithm for the digital restoration of degraded images in the presence of additive noise, modeled as a Pth order Markov random field, has been developed. This algorithm is superior to that of Hunt in two significant ways. First, the restriction of whiteness of the additive ~MMMMMMMMMHlllIllll ~MMMMMMMMMIIIIIIllll MMMMMMM~MMIIIIlIlIII ,MMMMMMM~MMIIIlIIIIII MMMMMMMMMMIIIIIIIIII MMMMMMMMMMHIIIIIIII MMMMMMMVMMIIIIIIIIII MMMMMMMMMMHIIIIIIII ]MMMMMMMMMMHIIIIIIII IMMMMMMMMMMH||||I||| lIIIIlOllIIIlOgOOIII IIIIHIIIIHIIIIIIII

MMMMMMZOOO||||||OZ÷- -*~ MMMZZMZOOOOI||OiOZ+ --4 MMMVMMZOOOHO|||OZ+- -+4 MMMVMMEEIOINllO||OZ+- - - 4 MMMMMMM|OOIIIgg||OM+- -+4 M M M M Z M Z I I B H | | | | | Z + - -+4 Z M Z Z Z Z Z e | O I I N | | I Z + - -+~ ZO|IZLOOO|HIIIBOZ+++M O|eiOOOO|OaOO|O|$Z+--+M2

IIIIHIIIIHIHIlIII IIIIlIIIIIHIIIIIlll

IIIIHIIHHIIIIOZ

IlIIllllllHlllllllI Illlllllllllllllllll IIIIHIIIIIIINIIlll

IIIlIlIIlllIIIIlIlIH IIIlHIIIIHIIIIllll IIIIHIIIIHIIlIIlII

|eeeeeoeHH|e||ez+--+za D|||ItI|||It||||OZM-IIIIllll IIIIIRI III1~

+ Zl Z + - - MZe

M÷- + Ze

IIIIHIIIIIDIIOIIIZM -MZl IIIIHIIIlIIIIIIOZ+--*Zi IIIIHIIIIIIIIIIgZ+--÷M0 IIeOIO~IOOO00~OeZM+---MI ZMZZMZZZMZZZZZZMM+-+M +÷+M÷M÷÷++M++++÷+--

....... ÷÷÷÷M÷÷M+---~÷÷÷+++MZZZZZZM+÷ ++++M÷+ZZOO|OaD~M+-

Fig. I. Due image.

Fig. 2. Distortedimage

---+

~ N ~m~

~

m

~

. m . . . ~ L~

)

~ m ~

~

o~mmm o.mmm (--mml +emmm

N

~ ~

m m m m

#

~

~

~ +

~

m m m m

~

m m m m

m m m m

m m l o

m m m m

m m m m m o m m m m m m -

#

I l l l l , + + _ + l l

I l l

5"

74

B. P. AGRAWALand V. K, JAIN

noise, crucial for Hunt's algorithm, has been removed. Second, the present algorithm requires considerably less core memory, of the order of N 3 rather than N 4. However, it is inferior to Hunt's algorithm in that it is relatively slow. The algorithm can be speeded up by using the fast Fourier transform (FFT) algorithm to compute the inverse of the matrix C'C, a block-circulam matrix and the product HF, a circular convolution. A very dramatic saving in computations can be realized in applications where a slight inaccuracy can be tolerated. This will be accomplished by exploiting the similarity of the block-Toeplitz matrices to the block-circulant matrices[16] to obtain a numerical inverse of the matrix (H'SH+Q) by the use of FFT methods. The implementation and testing of these ideas remains to be carried out, REFERENCES 1. B.R. Hunt, The application of constrained least squares estimation to image restoration by digital computer. IEEE Trans. Computers C22, 805-812 (Sept, 1973), 2. M. M. Sondhi, Image restoration: the remowd of spatially invariant degradations. Proc. IEEE, 60t7), 842-53 IJuly 1972). ~,. D. P. MacAdam, Digital image restoration by constrained deconvolution./. Opt. Soc. Am. 60, 1617-1627 (Dec. 1970). 4. C. L. Rino, Bandlimited image restoration by linear mean square estimation. J. Opt. Sot. Am. F, 547-553 (May 1969). 5. C. W. Helstorm, Image restoration by the method of least squares. Z Opt. Soc. Am. 57, 297-303, IMar. 19671. 6. J. L. Harris, Image evaluation and restoration. J. Opt. Soc. Am. 56, 569-574 (May 1966}. 7. J. W. Woods, Two-dimensional discrete Markovian fields. IEEE Trans. on In[o. Th., IT-18(2), 232-240 (Mar. 1972), 8. D. L. Phillips, A technique for the numerical solution of certain integral equalions of the first kind, J. Ass. ('otnpllt. Math. 9, 97-101 (1962L 9. S. Twomey, On the numerical solution of fredholm integral equations of the first kind by the inversion of linear system produced by quadrature. J. Assoc. Comp. Math. 10, 97-101 (Jan. t963). II). S. Twomey, ]'he application of numerical filtering to the solution of integral equatiuns encounter in indirect sensing measurements. J. Franklin Institute 279. 95-109 (Feb. 1965). 11. 13. R. Hunt, A theorem on the diificulty of numerical deconvolution, IEEE Trans. ,4udio Electr~;a(oustic~ AU 20, 94-05 tMar, 1972). 12. M. P. Ekstorm, A spectral characterization of ill-conditioning in numerical deconvolution. IEEE Trans, Audio Electracoustics AU 21, Ind. 4 (Aug. 1973). 13. 1,. M. Kutikov, The structure of matrices which are the inverse of the correlation matrices of random vector procc~e~ U.S.S.R. Comput. Math. and Phys. 25, 58-71 (1967). 14. H. Akaike, Block toeplitz matrix inversion. SIAM J. Appl. Math. 24, 234-241 (Mar. 1973). 15. B. P. Agrawal. Digital processing of two-dimensional signals limagesi. Ph.D. Dissertation, University of South Florida (I 974). 16. M. P. Ekstorm. An iterative-improvement approach to the numerical solution of ~ector toeptitz systems. IEEE Tran,s. Computers C-23, 320-325 (Mar. 1974). 17. (}. B. Anderson and T. S, Huang, Piecewise Yourier transformation for picture bandwidth compression. IEEE "I?ansuctions on Comm. l'ech. COM-I9{2), 133-141) (Apr. 1972).