State Space Approach for Two-Dimensional Data Smoothing

Copyright © IFAC Identification and System Parameter Estimation 1982. Washington D.C . . USA 1982

STATE SPACE APPROACH FOR TWO-DIMENSIONAL DATA SMOOTHING S. M. Hang and C. N. Shen Electrical, Computer, and Systems Engineering Department, Rensselaer Polytechnic Institute, Troy, New York 12181 , USA

Abstract. A state space approach for two-dimensional data smoothing using spline functions is developed in this paper. This approach not only keeps the results concise and simple, but it also provides the basis for obtaining a recursive algorithm. Minimization of a suitable criterion function results in an optimal solution which can be computed recursively if certain conditions are met. In case a recursive algorithm is feasible, certain spline functions can be used to represent the finite state recursive system. Conditions under which a particular spline function will yield a linear state model are derived. Specific examples are shown for the derivation of recursive smoothing algorithms using spline functions with linear state space models. I.

INTRODUCTION

11.

PROBLEM STATEMENT

A number of different approaches for 2-D data processing have been reported in the literature recently. These generally assume an underlying image model for the observed data for solving the prediction, filtering or smoothing problem [1-3).

The 2-D smoothing spline problem is defined as follows. Given a set of noise corrupted 2-D data Z = {z .. ; i = 1,2,---M, j = 1,2,---N} it is required f~ estimate the original function values and its derivatives from these discrete measurements.

In case a suitable two-dimensional model cannot be identified for an image, the data processing problem can be solved by minimizing a suitable criterion function, [4-8). The optimal solution for the criterion function used by Kim, (8), is nonrecursive in general and involves certain arbitrary functions that should satisfy specific differen-

The measurement process is described by i j

1,

M

1,

N

(1)

where f(~i' nj) are the sample values if the original function is assumed to have mixed partial derivatives up t o a certain order. Vij is random noise with zero mean and finite variance R ..•

tiability requirements, e.g., [ ;4s 2 ]2 = O. d lJd A

~J

The approximating function selected here is the 2-D spline function which is defined in wide sense; a piecewise fun c tion of continuous derivatives up t o certain order with each piece defined only in one grid, Fig. 1.

Kim, (8), used bicubic spline function bases for these arbitrar y functions and derived a sub-optimal recursive algorithm. This pap e r uses the s ame general criterion function as in (8), however, a wider class of possible signals is explored. Since a recursive algorithm can be considered as a finite state machine with feedback to update itself, a linear state space model for the recursive algorithm can always be found. Utilizing the useful fact that certain spline functions with finite support have linear state space models, we derive smoothing algorithms using such functions. The sufficient conditions are also derived under which the optimum solution can be made recursive and a particular spline function will yield a linear state model.

The reason for using spline function is that the approximating error analysis has been well developed in the numerical mathematics [4,5). Another advantage of this approach is that under certain conditions recursive computation can be realized, which makes the real time operation possible. The entire spline fun c tion in domain R as dictated in Figure 1 would be s( ~ , n )= sij( ~ , n )

~ i ":' ~..:. ~ i+l' nj < n ~ nj + 1 i = 0, 1 , .. • M-I, j =0, 1 ... N-l

379

(2)

S. M. Hang and C. N. Shen

380

The objective function which we wish to minimize is

M J

=

f

N

l:

l:

P

i=l j=l

ij

l:

6~

f611

o

0

[d 42Sij2 (\l.A)]

2d\ld A

(4)

d\l dA

is the one satisfying the Euler's equation,

M-I N-l

+

From the principle of calculus of variation, the optimum solution of

l:

d~dn

i=l j=l

d

(3a)

8

S ..

1J

(\l,A)

o

(5)

where

T -1

(~i,nj)l Rij [zij-s(~i,nj») (3b)

P

and

ij

's are weighting parameters as

o < P ij .::.

1

with proper boundary conditions [9). If only polynominal solutions are of concern, then the solution satisfies Equation (5) is equivalent to the one which satisfies the following equation,

C(s) is the smoothness measure function usually defined as a function of the high order derivatives of s(Cn). The first term in Equation (3a) represents the noise effect at the grid points and the second term is a measure of smoothness of the approximating function within the grid. The selection of C(s) and s(~,n) can be independent, but the optimum s(~,n) which minimizes a special

f J C(s)

[\lA A \l 1]· U . . 1J

where U .. 6 1J

~

(6)

[U~. U~. U:. U?) Tare under1J

1J

1J

1J

determined parameters. The general solution of Equation (6) is

dpd A is often unique.

Notice that s(~,n) would act as an interpolation function when P . . approaches to zero, 1J

and on the other extreme, when p .. approaches 1J

to one s(~,n) would become a smooth plane. Thus, the value of weighting parameter P ..

1J

controls the trade-off between the measurement error and the smoothness of the approximating function. Ill.

THE STATE SPACE MODEL OF BICUBIC SPLINE

The most often used splines are polynominal splines because of their simplicity in computation and the closed form solution. In general, a 2-D spline defined in one grid is described by a partial differential equation. However, only feW of the solutions satisfying some restrictions could be represented by linear discrete state space model defined on the grid points. The following example illustrates the formulation of linear state space model from a specified spline. This formulation procedure can be extended to the more general cases. The bicubic polynominal spline is the polynominal which minimizes the following smoothness measure C(s) in Eq. (4). For convenience, new notations of variables are introduced,

C, 1

for

A 6 n - n., J

for

p 6 c=

~ -

=~

0 '::'\l'::'

6~,

0 < A< 6n

~i .::. ~ .::. ~ i+l n. < II .::. nj+l J 6~ ~ ~i+l - C1 6n ~ nj+l - n

where ~l' ~2 ' ~l' and ~ 2 are arbitrary functions.

If we define the "state" as the vector

~ ~ ~] T , then the exact linear [ s, d\l ' dA ' d\ldA discrete model of Equation (7) defined on the end points of one grid can be obtained: (See Appendix). X'+ 1 l , J' +1

~

= 0,1 :

m

= 0,1.

(8a)

The linear discrete state model for a general 2-D spline is often not available. Because of the restriction of the uniform grid, variations of S( \l , A) propagate along the lines parallel to the S or II axis only. Consequently, the arbitrary functions could be functions of only one variable, 11 or A, so that these can be cancelled. Linearization or approximation may be used to derive the linear discrete model for those which do not meet above requirements. Since the finite state model is one of the essential condtions for recursive algorithm,

381

State Space Approach for Two-dimensional Data Smoothing the bicubic spline is the one, if not the only one, of the members which can be implemented recursively. IV.

(12b)

THE FORMULATION OF VECTOR PROCESSING

l-D estimation theory is well developed, nevertheless, the corresponding results in 2-D are not always possible. This complication occurs partially due to the choice of the support defined for the 2-D filters [10-12]. An estimator with support defined by Equation (8) usually does not have the optimum solution either in the sense of least square or quadratic performance index. To avoid this fundamental restriction, a vector processing model is constructed as follows. Define the global state quantity to be the vector along the column of region R, X ] Mj

T

j = 0,1;2 ... N

for bicubic splines. However, the result of the above problem can be easily extended to the general quadratic form,

M~

Jo

IOl~j+l

=

+

IlO~j

IOO~j

+

+

LO~j'

[U

Let

(9)

~j+1

F X. + ru. -J = -J

la],

(lla) (llb)

II

°

V ] Mj

J =

~

FOl

°

L

L

j=O

-1 R. =J

(lld)

°0

·Fo e

(llf)

(13)

-J

-1

~) Po

can be

(~ -~)

[z. _HX.]TR~l [~j - H X. ] -J = -J J = -J

U~ Q~l U.

-J

~J

diag[P

(14)

-J

ft,~ fOt,ll

°

a4s . . ] 2 [-1:..L a/aA 2

u~. Q-l u . . 1J

ij

Oj

-1 -1 R R Oj Plj lj

-1 ~ = -1 -1 . d1ag [(l-P Oj )QOj (l-P lj )Qlj

PMj

~~]

1J

-1 (l-P Mj )QMj]

The first term in Equation (14) is added to reduce the effect of the assumed initial condition. V.

To simplify the calculation and deduce the recursive solution, only the quadratic smoothness measure is considered, Eq. (6) gives =

H]

(14b)

LO = diag[r. r ... r]

C(s)d~dA

...

(14a)

0

(lle)

i° °

T

where ~ is the initial guess, and

lO lO rlO = diag[F lO F ... F ]

t,~ ft,ll

diag[H H

= -J

(~ -

N-l +

FOl.

0

T

The objective function (Equation (3» expressed in the global state form,

(llc)

0

FOO

[0 V V. V lj 2j -J ~

j=l

0

°FOO ° °FOO . °

T

z. = H X. + V.

+

° °

X.. X.+1 . X . . +1] 1J 1 ,J 1,J

ZMj]

N

and

°FOl

ij

-J

l l ] - l rO

[I

X.. X. 1 . X. . +1] C- l 1J 1+ ,J 1,J

[0 Zlj Z2j z. -J ~

H

(ll)

where

L

ij

and the corresponding observation equation becomes

=

_ [01]-1 [110 + 1= [I

C(s) d~dA

H ~ [1 o 0 0]

(10)

which could be rewritten as

0 [U

The state equation along the column strip is

~j+l

ft,ll

dlJdA (12) (12a)

THE SMOOTHING ALGORITHMS

We now turn to the development of the estimation algorithms. The closed form recursive formulas for the problem described in the previous sections does not always exist. There are some significant special cases that provide simple recursive computations as discussed below.

s.

382

M. Hang and C. N. Shen

For the sake of forward real time recursive computation, the unknown coefficients Qj'S can only be the function of previous states X. 's and parameters U. 's, i
~

-

(~j

T

-1 ~j

l!x'!j-1 - l!IQj-1)

(17)

Jj

~,~ ••• ,Qj-1

minX [JO +

HF X (~j - == -J-1 - l!IQj_1)

J. J

0

mi~

zi

[J 1 + ... [min U. Jjl ••. llll -4) -J-1 (15)

Since !j = X, !j -1 + IQj -1 . Then taking the gradient of Equation (17) we obtain the expression,

where

zi ~ J £,

=

(15a)

{z £, : 1,::, £, .::. j} H X)T -1 ) ( ~£, _ = _ £, R£, (~£, - ~!£,

T

+ Q£, -l

-1

or U. -J-1

(15b)

£, = 1,2, •.. j

Q£'-l Q£'-l

(18)

T~l

=

rT HT R~l = = =J

=,]

(~J.

- 1_!. I_.!J. -1) (19)

where (15c)

(19a)

and Because the preceeding solution is solved by using data The no tation lzi emphasizes that the measuremen ts are up t the j th stage. Let us denote the optimum estimate of !i j as X. ., that is, X.I. is

which minimizes J

o

-1

IJ

1

the optimum estimate of !i using the

J

dat~

zi,

Qj-1

and !j-1 should be

denoted as U. 11. and X. 11 . ' respectivel y . -J- J -J- J To compute the Qj-2 Equation (19) is substituted in Equation (17), and J. is simplified as a function of !j-1 only . J

(~j - l!x'!j-1)

J. J

T

-1

~j

(~j - l!I !j-1)

( 20)

~1' ~2' · ··~j·

and

Ri +1 l j x. 1 1.J

-1+

The relationship between Xilj c omes from Equation (8),

=

FX·1 IJ· + r -1 u · J· l

I

.,

J J

mi~

First, we compute U. I. 's from jth stage backth 1 J word t o 0 stage, and then the smoothed ~Ij

th

~j + l!

I

~-1

T

T

I li

(20a)

stage, we are going to solve j

j-2

J. 1 J-

we must complete all mi~

Ui l j' s and XOl j 's via Equation (15). The whol e procedure of computation is dictated in Figure 2 .

estimate

=

At j_1

zi.

A

~j

(16)

where the index i l j indicates that the quantitie s are the optimum solutions using In order to obtain X.

where

's are calculated forward by

using Equation (16).

-j-2

(J . 1 + JJ. ) J-

mi~ j -2 [(~j -1

+

T -1

li! j -1) ~ -1 (~j -1- li!j -1)

~-2 ,gj:2 ~-2 + (~ - liX, !j-1)

T -1 ~ (~j - l!x'!j_1) (21)

At the jth stage, U. 11. is obtained by minJ-

J

The resultant U. 2 1. is J-

imizing J. with respect to U. 1 J

J-

(z . - HX .)T -J =- J

R~l (_zJ. - =H_X ) J J

U . 21 J. -]-

T-. 1 {r T HT R-. 1 ( HFX) J-1 = = =,]-1 ~j-1 - = =-j-2

+ IT x,T + UT

-j-1

1"1-

1

~j-1

U -j-1

J

l

!i.j1(Zj - l! X, X,! j-2)}

(22) where

State Space Approach for Two-dimensional

(22a) Observing Equation (19) and (22), it seems that there is no simple closed recursive formula to compute U. I' 'so -1. J However, if the state equation (16) satisfies certain conditions, the closed recursive formulas for U' . 's or X. I . 's are possible. Sub-1. IJ -1. J stituting Equation (22) to Equation (21), j J J. 1 becomes function of X. 2 only. -J-

Repeat th the above procedure backward down to 0

stage, we will have Jd whi ch is a function !o only .. Then X ' can be obt a ined by minimizol ing with res~e c t t o XO . After that, XII' 'X . , .. . are then comp uted fon~ard via 21 J 6) J . , Equation (1 . This nonre c urS1.ve comput1.ng pro cedure can be applied to the general state model on llhich no re s tri c ti ons a re imposed.

Woods, J. W., and C. H. Radewan, "Kalman Filtering in Two-Dimensions," IEEE Trans. Infor. Theory, Vol. IT-23, pp. 473-482, July 1977.

4.

Ahlberg, J.H., E. N. Nilson, and J. L. Walsh, The Theory of Splines and Their Application, NY: Academic Press, 1967.

5.

Schultz, M. H., Spline Analysis, NJ: Prentice-Hal1, 1973.

6.

Hou, H.S., and H. C. Andrews, "Least Squares Image Restoration Using Spline Basis Functions," IEEE Trans. Comput., Vol. 26, pp. 856-873.

7.

Hou, H.S., and H. C. and H. C. Andrews, "Cubic Splines for Image Interpolation and Digital Filtering," IEEE Trans. Acoustic., Speech, Signal Processing, Vol. ASSP-26, pp. 508-517, Dec. 1978.

8.

Kim, C.S. and C. N. Shen, "Design of a Recursive Vector Processor Using Polynominal Splines," Proc. of 19th CDC, pp. 363-367, 1980.

9.

Gelfand, I.M., and S. V. Fomin, Calculus of Variations, NJ: Prentice-Hall, 1963.

DISCUSSION AND CONCLUSION

The general properties of 2-D smoothing sp1ines problem have been dis cussed in this paper. The principles and restrictions on forming a recursive algorithm of 2-D vector processor using splines are explored. It is found that the recursive on-line structure for general splines is not always possible. One special case which results in forward on -line recursive algorithm was developed. State space notation was used through as it not only keeps the results simple, compact, and unique but also provides a basis for recursive processing. Another factor involved in the formulation of recursive computation is the proper selection of the objective function which will be taken up in the future. Also, the adjustment of the weighting parameter Pij for adaptive processing is not included in this paper. The development of adaptive algorithms and the extension to splines with nonregular grids and/or nonlinear state model will be presented elsewhere.

10.

Marzetta, T.L., "Two-Dimensional Linear Prediction: Autocorrelation Arrays, Minimum-Phase Prediction Error Filters, and Reflection Coefficient Arrays," IEEE Trans. Acoust., Speech and Signal Processing, Vol. ASSP-28, pp. 725-733, Dec. 1980.

11.

Strintzis, M.G., "Comments on 'TwoDimensional Bayesian Estimate of Images," Proc. IEEE, Vol. 64, No. 8, pp. 12551257, August 1976.

12.

Barry, P.E., R. Gram, and C. R. Waters, "Two Dimensional Filtering--A State Estimator Approach," Proc. of 15th CDC, pp. 613-613, 1976.

In particular, it is hoped that the expositions will shed some light on the general problems of the recursive algorithms for 2-D smoothing splines, which does not appear in I-D data processing.

grid

REFERENCES 1.

Jain, A.K., "A semicausal Model for Recursive Filtering of Two-Dimensional Images," IEEE Trans. Comput., Vol. C-26, pp. 343-350, April 1977.

M- l

:'i

2.

Habibi, A., "Two-Dimensional Bayesian Estimate of Images," Proc. IEEE Vol. 60, No. 7, pp. 878-883, July 1972.

383

Smoothing

3.

J6

VI.

Dat~

H--t I I II II I.

Fig. 1

I ! ("r ' , 'N')

_ - - '_ _.....;._ _......:

The region of 2-D splines, R

S. M. Hang and C. N. Shen

384 APPENDIX For Eq. (8) in the text, we have

F

01

F

OO

1

0

I':,

0

1

0

I':,

I':,f,

0

0

1

0

1

0

0

0

1

1

I':,f,

0

0

0

1

0

0

0

0

1

0

0

0

1

I':,f,

I':,

0

1

0

n

10 F

n

0 n

(A-i)

1':,f,l':,n I':,

n

0

0

1

I':,f,

0

0

0

1

(A-2)

and 1':,31':, 3

1':,21':, 3

1':,31':,2

1':,21':,2

~ 36

~ 12

:Ql

~ 4

1':,21':,3

f, n

r

1':,21':,2

f, n

I':, 1':,2 f, n

1T

-6-

-4-

-2-

1':,31':,2

1':,21':,2

1':,31':,

1':,21':,

(A-3)

'=

f, n

----u-

Fig. 2

I':, 1':,3 f, n

12

f, n

-4-

1':,21':,2

I':, 1':,2

~ 4

.::1.::::!l 2

The Computation Scheme

f, n

-6-

f, n

-2-

1':,21':,

f, n

-2-

1':,f,l':,n