DECOMPOSITION AND INVERSION OF CONVOLUTION

OPERATORS

By

ZOHRA Z. MANSEUR

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL

OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT

OF THE REQUIREMENTS FOR THE DEGREE OF

DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1990

ACKNOWLEDGEMENTS

I would like to thank my advisor Dr. Gerhard X. Ritter, first for coming up

with such a great concept as the Image Algebra, and second for allowing me to be

part of his growing team. His guidance is also greatly appreciated.

This dissertation would not have been possible without the constant help of

my other advisor, Dr. David C. Wilson. His inquisitive mind encouraged me to

always look further. I would like to thank him for being there, available, patient,

and helpful. The generous use of his office space and facilities greatly helped in the

completion of this work.

Drs. James E. Keesling, Li-Chien Shen, and Joseph N. Wilson deserve thanks

for their service as members of my committee. I would like to thank Dr. Keesling

also for his willingness to serve as one of my references.

The financial support I received from the Mathematics Department in the form

of a teaching assistanship, and from Dr. Ritter in the from of a research assistanship

is gracefully ackowleged. I take this opportunity to thank the secretaries of the

Mathematics Department for caring and helping, when needed. I wish to thank also

Randy Fischer for all the help on computers.

I am deeply grateful to my parents, for providing me with a first class educa-

tion, despite the many difficulties they encountered in their life. My mother, Faiza

Houasnia, and my father, Tahar Zbiri, have taught me that a woman's independence

11 -

can be achieved mainly through her education. Special thanks go to my parents-in-

law Dhia Si-Ahmed, and Omar Manseur for their support and trust, and for raising

such a fine son.

Very sweet thoughts go to my husband and best friend, Rachid, for his under-

standing, and his patience throughout our marriage. He believed in me, as much

as my parents did. His constant help with computers gave me the privilege to have

almost all my technical questions answered. My children, Mehdi and Maya, made

my student life harder sometimes, but they also made it all worthwhile.

Finally, I wish to thank all my friends. They made Gainesville a very enjoyable

city to live in.

- 111 -

TABLE OF CONTENTS

page

ACKNOWLEDGEMENTS ................................................... ii

LIST OF SYMBOLS ................. ...... ........ ............. vi

A B ST R A CT ....................................................... ........... viii

C H A PT ER S .................................. .................. .............

1 INTRODUCTION ......................................................... 1

2 IMAGE ALGEBRA: AN OVERVIEW ...................................... 7

2.1 Introduction .......................... ............. ... ... 7

2.2 Definitons and Background .......................................... 8

2.2.1 Operands of the Image Algebra ................................... 8

2.2.2 Operations of the Image Algebra .................. .............. 10

2.3 Image Algebra and Linear Algebra ....................................... 11

3 SHIFT INVARIANT OPERATORS ........................................ 14

3.1 Decomposition Methods for General Operators ........................ 15

3.2 A Characterization of Operators Decomposable into a Special Form ....19

3.3 Decomposition methods for Symmetric Operators ................... 24

3.4 The LU Factorization ................................................ 31

3.5 Grobner Bases ................. ... .................................... 33

3.5.1 Terminology .................................. ................ 34

3.5.2 Application to Systems of Equations and Algorithm ............. 36

4 VARIANT OPERATORS ................................................ 39

iv -

4.1 Decomposition of General Operators ............................... 39

4.2 Decomposition of Symmetric Operators .............................. 46

5 EXAM PLES ............................................................. 50

6 OPERATOR INVERSION .............................................. 57

6.1 Introduction ....................................................... 57

6.2 The von Neumann Template ...........................................58

6.2.1 Preliminaries ................................................. 58

6.2.2 The Main Theorem ......................................... .62

6.3 The Characteristic Polynomial Method ...............................67

6.4 Inverse of General Shift Invariant Operators .........................69

7 CONCLUSION AND SUGGESTIONS FOR FUTURE RESEARCH .........74

REFERENCES .............................................................. 76

BIOGRAPHICAL SKETCH .................................................. 79

LIST OF SYMBOLS

Z:

R:

Rn:

F:

E, 0, C:

X, Y:

X\ Y:

x, y, z:

a:X - F:

a, b:

a + b:

FX:

t:Y --, FX:

ty:

S(ty):

r, s, t:

(FX)Y:

LX:

F[yl, Y2, --, Yn]:

the set of integers

the set of real numbers

n-dimensional Euclidean vector space

an arbitrary value set

is an element, is not an element, is a subset of

point sets (in particular n x n arrays)

the set difference of X and Y

points

a is an F valued image on X

images

addition of image a and image b

the set of F valued images on X

t is an F valued template from Y to X

the template t at the point y

the support of the template t at the point y

templates

the set of F valued templates from Y to X

the set (FX)X

the ring of multivariate polynomials over F

- vi -

R[x,y]/(xm 1,yn- 1):

Mn(F):

0, 77, 0:

Wn:

):

cir(co, cl, ..., Cn):

Cmn:

Fn:

Frnn:

the quotient ring of polynomials modulo x"' 1 and yn 1

the algebra of n x n matrices over F

isomorphisms

exp(-2ri/n), where i = (-l)

generalized convolution

an n x n circulant matrix

the set of circulant matrices

the one-dimensional Fourier matrix of order n

the two-dimensional Fourier matrix of order m x n

- vii -

Abstract of Dissertation Presented to the Graduate School

of the University of Florida in Partial Fulfillment of the

Requirements for the Degree of Doctor of Philosophy

DECOMPOSITION AND INVERSION OF CONVOLUTION

OPERATORS

By

Zohra Z. Manseur

December 1990

Chairman: Dr. Gerhard X. Ritter

Co-Chairman: Dr. David C. Wilson

Major Department: Mathematics

The image algebra, an algebraic structure developed by G. X. Ritter et. al.

specifically to meet the needs of digital image processing, will provide the mathemat-

ical setting for this investigation. The results presented in this dissertation evolved

as a direct result of questions that arose during the development of this image al-

gebra. A network of processors for a parallel computer architecture is modeled as a

subset of Rn. An image is represented as a function defined on such a subset. The

viii -

subsets of concern, here, will be rectangular arrays with m rows and n columns. A

template will be represented by a function from points in a finite rectangular array

into images. Basic operands and operations of the image algebra, such as addition,

multiplication and convolution, are defined. Relationships between image algebra,

matrix (linear) algebra, and polynomial algebra are described. The two main ques-

tions addressed in this dissertation are the problems of template decomposition and

template inversion.

Methods are established for a two-dimensional shift invariant convolution op-

erator defined on a rectangular array to be decomposed into sums and products of

operators of smaller size. The main result of this type is that 5 x 5 shift invariant

operators can be written as the sum and product of at most five 3 x 3 shift invariant

operators. The second is that a 5 x 5 operator satisfying certain symmetry condi-

tions can be written as the sum and product of at most three 3 x 3 operators. The

decomposition methods presented here are suitable for a machine with a pipeline

architecture.

An additional focus of this research will center on the problem of extending

results valid for shift invariant operators to non-shift invariant operators. In par-

ticular, necessary and sufficient conditions are given for a variant template to be

separable. A 5 x 5 variant template can always be written as the sum and product

of at most seven 3 x 3 templates.

The mean filter with five nonzero weights all equal to one and whose configu-

ration (or shape) is in the form of a cross is shown to be invertible when defined on

an m x m array if and only if neither 5 nor 6 divide m.

- ix

CHAPTER 1

INTRODUCTION

This dissertation is the result of work done in conjunction with an investiga-

tion of the structure of the image algebra, an algebraic structure for use in image

processing. The need for a unifying mathematical theory for image processing algo-

rithms specification, and a more powerful basis for an algebraically-based, high level

programming language, led to the development of the image algebra. It is an algebra

that operates on images, and whose operators reflect the types of transformations

commonly used in digital image processing [2]. It is capable of expressing all types

of transformations and algorithms commonly used in image processing, in terms of

its operators [28,38].

G. Matheron [18] and J. Serra [32] initiated the use of image algebra in image

processing. It was known then as mathematical morphology. Sternberg [33,34]

also investigated it independently. This algebra was based on the operations of

Minkowski addition and subtraction of sets in Rn [24]. Mathematical morphological

is about the study of shapes and patterns. The operations of Minkowski are often

referred to as dilation and erosion operations. Theses morphological operations

techniques can be applied to many image processing and image analysis problems.

However, they cannot, with the exceptions of very few simple cases, be used as a

basis for a general purpose algebraically based language for digital image processing.

Their weakness resides in the fact that they cannot easily express many commonly

-1-

-2-

used transformations such as linear image to image transformations and algorithms

[19]. In response to these limitations, and in order to establish an algebraic structure

capable of expressing most of the common image processing algorithms, G. X. Ritter

at the University of Florida began the development of a more general image algebra

capable of expressing all linear and morphological transformations [23,25,26,28].

This image algebra, will provide the mathematical setting for most of the re-

sults presented in this dissertation. In Chapter 2 basic operands and operations

are defined, and relationships between image algebra linear algebra and polynomial

algebra are discussed. Isomorphisms imbedding the linear algebra, and the poly-

nomial algebra into the image algebra give rise to the observation that a linear (or

generalized) convolution is equivalent to matrix multiplication, and the observation

that operator inversion (or deconvolution) is equivalent to matrix inversion. In the

case that an operator is circulant, operator inversion is equivalent to polynomial

inversion. In the usual implementation of the isomorphism between templates de-

fined on rectangular arrays and matrices, a shift invariant operator corresponds to a

block Toeplitz matrix with Toeplitz blocks, and a circulant operator corresponds to

a block circulant matrix with circulant blocks. The image algebra allows us to give a

precise formulation of the decomposition methods for non-shift invariant operators.

Recall that some frequently used non-shift-invariant operators include the Fourier

Transform or the Gabor Transform [7,11].

The use of parallel processing has come to play an increasingly important role

in image processing during the past few years. Efficient implementation of linear

convolution is an important aspect of any hardware environment. Some classical ex-

amples of two-dimensional operators include the Gaussian mean filter, the Laplacian

filter, and the Discrete Fourier Transform.

-3-

Template decomposition plays a fundamental role in image processing algo-

rithm optimization, since it provides means for reducing the cost in computation

time, and therefore increases the comptutational efficiency of image processing al-

gorithms. This goal can be achieved in two ways, either by reducing the number

of arithmetic computations in an algorithm or by restructuring the algorithms to

match the structure of a special image processing architecture. One of the goals of

this work is to present methods for the decomposition of a two-dimensional shift-

invariant convolution operator in a form suitable for a parallel image processing

machine with hardware similar to the PIPE developed at the Center for Manufac-

turing Engineering of the National Bureau of Standards. As mentioned by O'Leary

[22], a machine of this type is able to efficiently apply a 3 x 3 convolution operator

to a 256 x 256 image, but has a more difficult time applying a 5 x 5 operator.

Another example of a pipeline computer is ERIM's CytoComputer [34] which can

deal only with templates of size 3 x 3 or smaller on each pipeline stage. Thus, un-

less Fourier Transform methods are to be used, operators of size larger than 3 x 3

must be decomposed into sums and products of 3 x 3 operators to fit the hardware

restriction. O'Leary showed that if an operator (viewed as a matrix) has rank n,

then it can be written as the sum of at most n rank 1 (i.e. separable) operators.

An immediate consequence of the rank one method described by O'Leary in [22] is

that any arbitrary 5 x 5 shift-invariant operator can be written in the form

t = (tl E t2) + (t3 E t4) + (t5 E t6) + (t7 E t8) + (t9 g t0i),

where each ti is a 3 x 3 shift-invariant operator and e denotes the convolution

operation.

While the results on shift invariant operators are presented in Chapter 3 as

theorems concerning the decomposition of polynomials of two variables into sums

and products of lower degree polynomials, the theorems discussed in Chapter 4 will

-4-

all be phrased in terms of sums and E operations on templates. In Chapter 3, it is

proved that an arbitrary 5 x 5 shift-invariant template t defined on Z x Z can be

decomposed into the form

t = (tl t2) + (t3 t4) +t5,

where each ti is a 3 x 3 shift-invariant template. Thus, the number of 3 x 3

templates needed to decompose t is actually half of that indicated by the rank

method. Example 5 of Chapter 5 shows that this decomposition is the best possible,

in the sense that it may not always be possible to decompose an operator into the

sum and product of fewer than five 3 x 3 templates.

The most commonly used template in the decompositions is the 3 x 3 square

template, better known as the Moore neighborhood. It has the following form:

Since most templates used in image processing are symmetric or skew symmet-

ric with respect to the x or the y axis, a discussion of the problem of decomposing

convolution operators exhibiting some type of symmetry is included in Chapter 3

and Chapter 4.

An second approach to the decomposition of two dimensional operators based

on the lower-upper triangular (LU) factorization of a rectangular matrix, will result

into the finite sum of separable operators. A comparison of the polynomial method

and the LU factorization methods is given in Chapter 3 Section 4.

-5-

Results concerning the decomposition of arbitrary operators, (i.e. those which

may not be shift invariant) are presented in Chapter 4. Necessary and sufficient

conditions for a variant operator to be separable, i.e. splittable into the convolution

of a horizontal operator and a vertical operator, are established. Since not all

templates are separable, we will show that every variant template can be written as

the sum of separable templates.

We conclude by pointing out that template decomposition is not only necessary

in the case where a special image processing device cannot handle large templates

but is also desirable when the problem of efficient decomposition is concerned in

real-time applications.

Inverse problems have come to play a central role in modern applied mathemat-

ics. While classical linear and non-linear inverse problems abound in mathematical

physics [31], they are also important in such imaging areas as tomography [31],

remote sensing [31], and restoration [30]. While Rosenfeld and Kak [30] showed

how the Wiener filter and least square methods can be used to restore an image,

they also explained the equivalence between these techniques and the problem of

inverting a block circulant matrix with circulant blocks. In 1973, G. E. Trapp [36]

showed how the Discrete Fourier Transform could be used to diagonalize and invert

a matrix that was either circulant or block circulant with circulant blocks. While the

algebraic relationship between circulant matrices and polynomials was completely

formulated by J. P. Davis [6], it was P. D. Gader and G. X .Ritter [9,27] who made

the connection between polynomials and circulant templates. In particular, it was

found that a circulant template defined on a rectangular array with m rows and

n columns is invertible if and only if its corresponding polynomial p(x, y) has the

property that p(wj,w ) 7 0, for all 0 < j < n, and 0 < k < m. (The symbol wn

denotes the root of the unity, exp(2ri/n).)

-6-

An application of this method is that the usual 3 x 3 mean filter defined on a

rectangular array with m rows and n columns is invertible if and only if the number 3

does not divide either m or n. Thus for example, the mean filter defined on an array

with 512 rows and 512 columns will be invertible, while the mean filter defined on

an array with 512 rows and 240 columns will not be invertible. Gader [8] asked if a

similar result can be obtained for the von Neumann mean filter. The von Neumann

mean filter is the template defined by

1

t = 1 <1> 1

1

The main result of Chapter 6 is to prove that if t is defined on a square array

with m rows and m columns, then t is invertible if and only if neither 5 nor 6 divide

m.

A variety of examples are presented in chapter 4 to both illustrate the utility

of the methods developed and to demonstrate the fact that some theorems are best

possible.

CHAPTER 2

IMAGE ALGEBRA: AN OVERVIEW

2.1. Introduction

The following chapter contains a brief review of the fundamental concepts and

notation of the image algebra. We will need these ideas to extend results valid for

shift-invariant operators to more general operators. For the purpose of this disser-

tation, only basic definitions and operations are given. The relationship between

image algebra and linear algebra will also be described in this chapter. For a full

and detailed description of all image algebra operands and operations, see [29].

The image algebra is a heterogeneous algebraic structure specially designed to

meet the needs of digital image processing. It is used as a model and tool for the

development of local, parallel algorithms. Several commonly used image processing

transformations, such as general convolutions, Discrete Fourier Transform, and edge

detection techniques can be expressed in terms of the image algebra.

The use of image algebra in digital image processing was initiated by Serra [32],

Miller [19], and Sternberg [33]. It was Ritter [24] who showed that their algebras

were all equivalent and based on the morphological operations of Minkowski.

An image algebra is an algebra whose operands are images and subimages

(or neighborhoods). It deals with six basic type of operands, namely, value sets,

-7-

-8-

coordinate sets (or point sets), the elements of the value and the coordinate sets,

images and templates.

2.2. Definitions and Background

2.2.1. Operands of the Image Algebra

A value set is the set of values that an image can take. It can be any semi-group

with a zero element. The most commonly used ones in image processing are the set

of positive integers Z+, integers Z, rational numbers Q, real numbers R, positive

real numbers R+, or complex numbers C. An unspecified value set will be denoted

by F.

A coordinate set (or a point set) is a subset of an n-dimensional Euclidean space,

Rn, for some n. It is commonly denoted by X or Y, and the elements of such sets

are denoted by lower case letters. The familiar point sets are the rectangular and

hexagonal arrays.

DEFINITION 2.1. Let X, and F be a point and a value set, respectively. An F valued

image a on X is a function a : X F.

Thus, the graph of an F valued image a on X is of the form

a = {(x,a(x)): a(x) e F, for all x e X}.

The set X is called the set of image points of a, and the range of the function a

is called the set of image values of a. The pair (x, a(x)) is called a picture element

or a pixel x the pixel location and a(x) the pizel (or gray) value We will denote

-9-

the set of all F valued images on X by FX. We make no distinction between an

image and its graph.

DEFINITION 2.2. An image a : X -- F has finite support on X if a(x) = 0 for only

a finite number of elements x E X.

Another basic, but very powerful tool of the image algebra, is the generalized

template.

DEFINITION 2.3. Let X and Y be two coordinate sets, and let F be a value set. A

generalized F valued template t from Y to X is a function t : Y FX

Thus, for each y E Y, t(y) e FX, or equivalently, t(y) is an F valued image

on X. For notational convenience, we define ty t(y). Thus,

ty = {(x,ty(x)) :x e X}.

The sets Y and X are called the target domain and range space of t, respectively.

The point y is called the target (or domain) point of the template t, and the values

ty(x) are called the weights of the template t at y. Note that the set of all F valued

templates from Y to X can be denoted by (FX).

If t is a template from Y to X, then the set

S(ty) = {x E X : ty(x) 0}

is called the support of ty.

If t is an F valued template from X to X, and X is a subset of Rn, then

t is called translation invariant (or shift-invariant) if and only if for each triple

x, y, z E X, with x + z and y + z e X, we have that

10 -

ty(x) = ty+z(x + z).

x

Note that a translation invariant template must be an element of (FX) Invariant

operators on Z x Z are commonly expressed in terms of polynomials of two vari-

ables. For brevity, we will frequently refer to shift-invariant templates as invariant

templates. A template which is not necessarily translation invariant is called trans-

lation variant or, simply, a variant template. Translation invariant templates occur

naturally in digital image processing.

2.2.2. Operations of the Image Algebra

The basic operations on and between F valued images are naturally derived

from the algebraic structure of the value set F.

Let X be a subset of Rn. Suppose a E RX and t E (RX)X.

Addition on images is defined as follows: If a, b E FX, then

a + b {(x, c(x)) : c(x) = a(x) + b(x), x X}.

Higher level operations are the ones that involve operations between templates

and images, and between templates only.

The addition of two templates is defined pointwise. If s and t E (RX)X, then

we have

(s + t)y() = Sy(x) + ty(x).

DEFINITION 2.4. The generalized convolution of an image a with a template t is

defined by

at {(y, b(y)): b(y)= 3 a(x)ty(x), y X}.

xeX

11 -

Linear convolution plays a fundamental role in image processing. It is involved

in as many important examples as the Discrete Fourier Transform, the Laplacian,

the mean or average filter and the Gaussian mean filter.

DEFINITION 2.5. If s and t are templates on X, then we define the generalized

convolution of the two templates as the template r = s E t by defining each image

function ry by the rule

ry = {(z,ry(z)) : ry(z)= 7 ty(x)sx(z), where z e X}.

xeX

Note that r can be viewed as a generalization of the usual notion of the com-

position of two convolution operators. If the templates r and s are translation

invariant, then except for values near the boundary, the previous definition agrees

with the usual definition of polynomial product.

Note also that if s and t are two invariant templates, then r would be an

invariant template too. Defining r at an arbitrary y E Y is sufficient to define the

template everywhere.

Many other image operations are described in detail in Ritter et al [29]. A

precise investigation of the linear operator 6 can also be found in Gader [8], and an

extensive study of other non-linear template operations can be found in Davidson

[5] and Li [16].

2.3. Image Algebra and Linear Algebra

If X is a finite rectangular subset of the plane with m rows and n columns,

then it can be linearly ordered left to right and row by row. Thus, we can write X =

12 -

{x1,X2,...,Xmn}. Let LX denote the set of templates on X. (i.e.

LX = (FX)X.)

Let (Mmn,,+,*) denote the ring of mn x mn matrices with entries from F

under matrix addition and multiplication. For any template t, we define a matrix

Mt = (mij) where mij = txj(x). For the sake of notational convenience we will

write tji for txj(xi).

Define the mapping : LX -- Mmn by 0(t) = Mt.

The next Theorem was proved by Ritter and Gader [8]. It shows that there is

an embedding of linear algebra into image algebra.

THEOREM 2.6. The mapping ) is an isomorphism from the ring (LX, +, e) onto

the ring (Mmn,+,*). That is, if s,t E LX, then

1.) 0(s + t) = 0(s) + 0(t) or Ms+t = Ms + Mt,

2.) 0(s e t) = 0(s)O(t) or MsE t = MsMt,

3.) b is one-to-one and onto.

This theorem clearly states that template inversion or deconvolution is equiv-

alent to matrix inversion. Actually, a more powerful implication of this theorem is

that any tool available in linear algebra can be directly applicable to problems in

image algebra.

DEFINITION 2.7. Let X be an m x n coordinate set. We say that the mapping :

X -> X is a circulant translation if and only if 4 is of the form

(x + h) = (x + h)(mod(m, n)), for some h E X.

13 -

DEFINITION 2.8. We say that t E (FX)X is circulant if and only if for every cir-

culant translation q, the equation tx(y) = t(x)(O(y)) holds.

This last definition shows that a circulant template is completely determined

if it is defined at only one point. Translation invariant and circulant templates are

used in the implementation of the convolution and therefore, are directly related to

the Discrete Fourier Transform.

For the purposes of this dissertation, X will usually be a subset of Z x Z,

where Z denotes the set of integers. All the templates we will consider will have

finite support contained in a set of the form

Xm,n = {(i,j) : |il < m and Ijj < n}.

Therefore, a (2m+ 1) x (2n+ 1) template on Z x Z will have support contained

in Xm,n and will be displayed as a matrix of the form

t-m,-n ... t-m,O ... t-m,n

t = tO,-n ... < tO,0 > ... tO,n

tm.-n ... tm,0 ... tm,n

If the template t is defined as above, then it is understood that all weights

tij = 0, for all (i,j) such that |il > m or Ijl > n. The element < too > will denote

the weight at the center pixel location (i.e. tx(x) = to0,). Note that the indexing

here differs by a shift from the one given above.

CHAPTER 3

SHIFT INVARIANT OPERATORS

The use of parallel processing has been increasing for the past years. Linear

convolution is widely used in image processing. It consists of applying an operation

between a template and a given input image, pixelwise, to yield an output image.

One of the main reasons for template decomposition is that some current image

processors can only directly implement operations on very small templates Most

transforms are not able to be implemented directly on a parallel architecture. Con-

sider, for example, a parallel image processing machine with hardware similar to

the PIPE, developed at the Center for Manufacturing Engineering of the National

Bureau of Standards. This machine is capable of evaluating a small template like

a 3 x 3,but cannot easily evaluate a template of larger size.

The problem of template decomposition is to try to find a sequence of tem-

plates t1,t2,...,tn smaller in size than the given template t, such that when we

apply an operation on template t and image a, it will be equivalent to sequentially

applying operations between t1,t2, ..., t to the image. For example, in the case of

the neighborhood array processor, it will consist of decomposing the template into a

sequence of factors, where each factor is directly implementable on the architecture.

The motivation behind all the work on template decomposition is to increase

the speed and reduce the convolution computation time. If a convolution operator

is large in size, then the cost in computation time can be prohibitive. This compu-

14 -

15-

tational complexity can be reduced if the given template in decomposed into sums

and convolutions of smaller size templates.

The following discussion will clarify the use of template decomposition. Sup-

pose a template t has the following decomposition:

t= (i=hri) + ( -qsj).

The convolution of an image a together with the template t will give

aet=a (=ri + sj=lsj)

= [(...((ae rl) e r2)...) rh] + [(...((a sl)E s2)...) sk]

=b e c,

where b = [(...((a rl) E r2)...) rh], and c = [(...((aE si) E s2)...) E sk].

When a t is computed, the image b is computed first, the result stored, then

the image c is computed and the sum of the two is taken.

3.1. Decomposition Methods for General Operators

Templates come in different weights, size or shape, depending on the image

processing application. A special classes of invariant templates are the two dimen-

sional rectangular templates. The support of a rectangular template is a rectangular

array.

16 -

While a 5 x 5 operator is usually defined as a matrix of the form

t-2,2

t-2,1

t-2,0

t-2,-l

t-2,-2

t-1,2

t-1,0

t-1,-1

t-1,-2

t02

to1

< 00 >

t0,-1

t0,-2

t12

tll

tio

tl,-1

tl,-2

t22

t21

t20

2,-1

t2,-2

we will express all the theorems, corollaries, and propositions in terms of polynomi-

als. To avoid negative exponents, we begin all indexing with the integer zero. This

implies, if n = m = 4, that the remainder term r(x, y) is shifted, so that its value

at the center pixel location will be represented by r22. Thus, in the decomposition

methods presented in Sections 2, 3, and 4, it will never be necessary to shift the

image before applying the operator r(x,y).

THEOREM 3.1. If t(x,y) is a polynomial defined by the formula

m n

t(x,y)= E Etijx y

i=0 j=0

then there exist five polynomials pl(x, y), p2(x, y), ql(x, y), q2(x, y) and r(x,y) such

that

t(x,y) = Pl(x,y)ql(x,y) + p2(x,y)q2(x,y) + r(x,y),

2 2

Pi (X)yY) = Epi jjxiyj,

i=o j=o

2 2

p2(x, y) = ET P2ijxy

?=o j=0

?n-2 n-2

q, (XI Y)= T T qlijx i yi

i=0 j=0

m-2 n-2

q2(X,y) = C q2ij xyx7

i=0 j=0

where

-17-

and

m-1 n-1

r(x, y)= Z : rijXy3.

i=1 j=1

PROOF:

If a(y) = E 00 ti, then by the Fundamental Theorem of Algebra, there are

two polynomials al(y), and a2(y) with deg(al(y)) < 2 and deg(a2(y)) < n 2,

such that a(y) = al(y)a2(y). Similarly, if c(y) = E niY, then there are two

polynomials cl(y), and c2(y) with deg(cl(y)) < 2 and deg(c2(y)) < n- 2, such that

c(y) = cl(y)c2(y).

If we choose

pl(x,y) = al(y)+x2c (y),

qi(x,,y) = a2(y) + m-2c2(y),

then the polynomial

u(x,y)= t(x,y) pl(x, y)ql(x,y)

m n

= E Uijx

i=0 j=0

will have the property that if uij is nonzero, then 1 < i < m 1.

If b(x) = E-i1 uioxi, then the polynomial b(x) can be decomposed as xbl(x),

where deg(bl (x) < m 2. Similarly, if d(x) = E n-T1 2z then d(x) can be

decomposed as xdl(x), where deg(di(x)) < m 2.

Let

P2(x,y) = x + xy2, and

q2(x, ) = b(x) + dl(x)yn-2.

If

r(x, y) = u((, y) p2(x, y)q2(x, y)

?m n

= 0 rijx ,

i=O i=O

18 -

then when rij is nonzero, i and j will be such that 1 < i < m 1, and 1 5 j < n- 3.

Q.E.D.

COROLLARY 3.2. If t(x,y) is a polynomial defined by the formula

4 4

t(, y) = E tijiy j,

i=0 j=0

then there exist five polynomials pl(x,y), p2(x,y), ql(x,y), q2(x, y), and r(x,y) such

that

t(x, y) = p1(x, y)ql(x, y) + P2(x, y)q2(x, y) + r(x, y),

where

2 2 2 2

Pl(x,y) = E EPlijxiyj' q(x,y) = E E qlijxiyj

i=o j=o i=0 j=

2 2 2 2

P2(x,y) = E EP2ijxi y, q2(x y) = q2iji,

i=oj=0 i=0 j=0

and

3 3

r(x,y) = rij2x2yJ.

i=1 j=1

Remark 3.3. An immediate consequence of Theorem 3.1 is that if t(x, y) is a poly-

nomial defined by t(x, y) = E=o j o tijgx2y3, then there exist n 1 polynomials

tl(x,y), t2(, y),..., tn-2(, y), and r(, y) such that

t(x, y) = tl(x, y) + ... + tn-2(x, y) + r(x,y),

where each ti(x, y) is a polynomial which can be written as a product of polynomials

of the form

2 2

p(x, ) = E pijxij,

i=0 j=0

19 -

and

3 3

r(x,y)= Z j rijxyJ.

i=1 j=1

3.2. A Characterization of Operators Decomposable into a Special Form

THEOREM 3.4. Let t(x,y) be the polynomial defined by :

m n

t (xy) = tij xiy ,

i=0 j=0

where too, ton, tmo and tmn are all nonzero.

Let

n mn

a(y)= toiy Y b(x)= tinxi

i=0 i=0

n m

c(y) = tmi,,y, and d((x) = tioxi.

i=0 i=0

Let al,a2,...,an be the roots of a(y); /1, 2, .., m be the roots of b(x);

71,72, ...,7n be the roots of c(y); and 61,62,..., 6m be the roots of d(x).

If there exists an ordering of the roots so that we have both equations

ala26162 = 01027172, and

a3...Cn.. ...S-6m = 3- ..m73-..7n

(or both equations

ala2012 = 71Y251S2, and

a3 ... an3 ... m = 73---7Yn3...mn),

satisfied, then there exist three polynomials p(x,y), q(x,y) and r(x,y) such that

t(x, y) = p(x, y)q(x, y)+ r(x, y),

- 20 -

where

2 2

p(x, y) = P: Pj y,

i=0 j=O

m-2n-2

q(x, y)= qijjxy3,

i=0 j=0

and

m-1 n-1

r(x,y) = Z rijxyj.

i=1 j=1

Conversely, if

t(x, y) = p(x, y)q(x, y) + r(x, y),

where p(x,y),q(x,y) and r(x,y) are defined as above, then there exists an ordering

of the roots so that we have both equations

ala26162 = /1 27Y172, and

a3...an3.---.m &= 3... ,-73---7n

(or both equations

Qala212 = 71726162, and

a3...an03--- ... = 73---7n53.---m)

satisfied.

- 21 -

PROOF:

Factor a(y), b(x),c(y) and d(x) to obtain

a(y) = ton(y + ... + (tol/ton)Y + (too/tOn))

= ton(y2 (al + a2)y + ala2)(yn-2 + .. + (-)na3...an), (1)

b(x) = tmn(xm + ... + (tln/tmn)x + (ton/tmn))

= (x2 (/1 + 32)x + 0132)(tmnxm-2 + + (-l)mtmn33...3m), (2)

c(y) = imn(y + .- + (tml/tmn)y + (tmo/tmn))

= (y2 (71 + 72)Y + 7172)(tninnn-2 + + (-1)ntmn73---7n) (3)

d(x) = tmo(xm + ... + (tlo/tmo) + (too/tmo))

= tmO(X2 (61 + 62)x + 6152)(xm-2 ...- + (-l)m 3...S), (4)

then we have the relations :

(-1)nalQ2...an = too/tOn, (1')

(-1)m1/32...Am = tOn/tmn, (2')

(-1)n7172..-n = tmo/tmn, and (3')

(-i)m162 ...Sm = t00/tno0- (4')

If we let

al(y) = P1P2yy2 /1/2('1 + 2)Y + 12 1\2,

and

a2(y) = (-1)mtmn3....my n-2 + ... + (-)m+ntmn 3...ma33...can,

then by substituting formula (2') in equation (1), it will be true that

a(y) = al(y)a2(y).

-22 -

If

dl(x) = '71722 Y7172(61 + 62)x + 71726162,

and

d2(y) = (-1)ntmny3...-^nXm-2 + ... + (-)m+ntlmn3"3..6n3...Sm,

then by substituting formula (3') into equation (4), it will be true that

d(x) = dl(x)d2(x).

Similarly, we can rewrite formulas (2) and (3) as

b(x) = b(x)b2(x) and c(y) = c(y)c2(y).

To define the polynomial

p(x, y) = EO Ej=O Pijxiy3

let p(l,1) = 0, poi be the coefficient of y2 in al(y), p2i be the coefficient of y2 in

cl(y), Pio be the coefficient of x2 in d1(x), and Pi2 be the coefficient of xi in bi(x).

To define the polynomial

q(x,y) = m-2 n-2 ij yj,

q y= i=0 j=0 qijx y ,

let qoj be the coefficient of y3 in a2(y), qm-2,j be the coefficient of yJ in c2(y), qio

be the coefficient of x' in d2(x), and qi,n-2 be the coefficient of x' in b2(x), and set

all other coefficients equal to zero.

If we set

r(x, y) = t(x, y) p(x, y)q(x, y)

7-, n

i=0 j=0

then, it follows from the conditions on the roots that if rij is nonzero, then 1 < i <

- 23

m-l1and 1

The converse can be proved by observing that if

t(x, y) = p(x, y)q(x, y) + r(x, y),

where

r(x,y) = E n-1 rij 2y',

then equations (1)-(4) can again be used to show that either

crl2S162 = 127Y17Y2, and

3...On63...6mn = 333...3mY3-..7n

or

alQc2/1 2 = 71726162 and

a3...nC3...flmm = 73...-7n3-...m.

Q.E.D.

COROLLARY 3.5. If t(x, y) is a polynomial defined by the formula

4 4

t(x, y)= tijx yj,

i=0 j=0

and t(x, y) satifies the condition of Theorem 3.4, then there exist three polynomials

p(x,y), q(x,y), and r(x,y) such that

t(x, y) = p(x, y)q(x, y) + r(x, y),

wh ere

2 2 2 2

p(x,y) = ~ Pijx'y, q(x, y) = E qijx2yJ,

i=0j=0 i=oj=0

- 24 -

and

3 3

r(x,y) = E rijxiy.

i=1 j=1

3.3. Decomposition Methods for Symmetric Operators

Most operators used in image processing are symmetric or skew-symmetric

with respect to the x or y axis. The following section deals primarily with operators

exhibiting some kind of symmetry.

DEFINITION 3.6. A polynomial p(x) = anxn + anl-n-1 + ... + alx + ag, is said

to satisfy the symmetric property with respect to n, if ai = ani for all i = 0,...,n.

Remark 3.7. We do not necessarily assume that the degree of p(x) is equal to n.

For example, if p(x) = x, then p(x) satisfies the symmetric property with respect to

the integer 2 .

DEFINITION 3.8. A polynomial p(x) is said to satisfy the skew symmetric property

with respect to n, if ai = -ani for all i = 0, ..., n.

Remark 3.9. If p(x) satisfies the symmetric property and has even degree, then the

exponents of xi, and xn-i are of the same parity. Therefore, if x = 0 is a root of

multiplicity k, then p(x) = xkpl(x), where pi(x) is of even degree.

25 -

PROPOSITION 3.10. If p(x) is a polynomial satisfying the symmetric, or the skew

symmetric property with respect to some integer n, then a non-zero number a is a

root of p(x) if and only if 1/a is a root of p(x).

PROOF:

Since p(1/a) = (l/oan)p(a) and a is a non zero root, the polynomial p(a) =

0 if and only if the polynomial p(l/a) = 0.

Q.E.D.

PROPOSITION 3.11. If p(x) is a polynomial of even degree n satisfying the symmet-

ric property, then there exist two numbers ki and k2 such that:

p(x) = axklql(x)q2(x)...q2(),

where a is the leading coefficient ofp(x), and each qj(x) is of the form x2 + bx + 1.

PROOF:

If a is the coefficient of the highest degree of x in p(x), and k1 is the largest

integer such that xki divides p(x), then p(x) = axklpl(x), where pl(x) satisfies the

symmetric property. Since pi(x) is symmetric, and of even degree, we know by

Proposition 3.10 that each root of pl(x) can be paired with its reciprocal. Thus,

pl(x) has n/2 factors of the form (x + c)(x + 1/a) = x2 + (a + 1/a)x + 1. Thus,

pl(x) = (x2 + blx + 1)...(x2 + b(n-k2)/2x + 1).

Q.E.D.

PROPOSITION 3.12. If p(x) = aixz is a polynomial satisfying the skew sym-

metric property with respect to an even integer n, then there exists a polynomial

p1(x) such that p(x) = (x2 l)pl(s), where pl(x) satisfies the symmetric property.

- 26 -

PROOF:

Since ai = -ani, where n is even,

2-n n/2

p(1) = Ei=oai = E oai E=oai = 0.

Similarly, p(-l) = 0. Thus, x2 1 is a factor of p(x).

We must now show that if p(x) = (X2 1)pl(x), then Pl(x) satisfies the sym-

metric property. If

P(x) = bn-2Xn-2 + ... + bix + bo,

then note that

ao = -bo, a = -bi, a- = bn-3, and an = bn-2-

Since p(x) satisfies the skew symmetric property, bn-2 = bo, and bn-3 = bl.

Note also that ak = bk-2 bk, for all k = 2,..., n 2.

We must now show that bk = bn-2-k. If we assume inductively that bk-2 =

bn-k, then, since p(x) satisfies the skew symmetric property,

bk-2 bk = ak = -an-k = -(bn-k-2 bn-k),

so that bk = bn-2-k. Thus, pl(x) is symmetric.

Q.E.D.

COROLLARY 3.13. If p(x) is a polynomial satisfying the skew symmetric property

with respect to an even integer n, then there exist two numbers k1 and k2 such that

p(x) = a(x2 J)xkzql(x)q2(x)...qk,(x),

where a is the leading coefficient of p(x), and qj(x) is of the form (x2 + bjx+ 1), for

j = 1,..., k2.

27 -

PROOF:

If zero is a root of p(x), then let k1 be its multiplicity. By Proposition 3.12 and

Remark 3.9, there is a symmetric polynomial pl(x), such that

p(x) = anxkl(x2 _- )p1(x).

Since n is even, and pl(x) satisfies the symmetric property with respect to n,

Proposition 3.11 implies that pl(x) factors into a product of quadratic polynomials

of the form x2 + bx + 1.

Q.E.D.

Remark 3.14. If

m n

p(x, ) = 7 Y tjx yj

i=0 j=0

then p(x, y) can be written in the form

=0o sX(Ej=o0 ijYj) = E-=o xipi(y).

The polynomial p(x, y) can also be written in the form

lj=o y( Zo tijx2) = lj=o YJqj(x).

DEFINITION 3.15. If p(x,y) =ETon0 j=ot = ix 0 oipi(y) then the poly-

nomial p(x, y) is said to satisfy the symmetric (respectively the skew symmetric)

property with respect to y, if pi(y) satisfies the symmetric (respectively the skew

symmetric) property, for all i = 0, ...,n.

A similar definition can be given for the variable x.

28 -

PROPOSITION 3.16. If p(x,y) is any polynomial, then

p(x,y) = pl(x,y) + p2(x,y),

where Pl (x, y) satisfies the symmetric property with respect to y, and p2(x, y) satisfies

the skew symmetric property with respect to y.

PROOF:

Let pl(y)= Zj=0 aj+an-j and 2(y) = j=0 ajy-an-

Q.E.D.

A similar proposition can be given for the variable x.

PROPOSITION 3.17. If t(x,y) = =0 tjx t y3 has the property that

t00 = ton = tmO = tmn = 0; and to0 = ktol, trn,n-1 = ktm-l,n, tm-1,0 = ktml,

and to,n-1 = ktin, where k = +1 or -1, then there exist three polynomials p(x,y),

q(x,y), and r(x,y) such that :

t(x, y) = p(x, y)q(x, y) + r(x, y),

where

2 2 m-2 n-2

p(x,y) = pijyxy, q(x,y)= q ijx2yJ,

i=0 j=O i=0 j=o

and

77m-1 n-1

r(xy) = i= ri j=

i=l j=l

- 29 -

PROOF:

If k = 1, let

p(x,y) = + y + xy2 + 2y,

q(x, y) = [t1 + ... + 0,n-yn-2] + [lnyn-2 + + tm-1,n m-2 n2

+ [tmlxm-2 + + m,n-2 n-3] + [t20x + -. + tm-2, -3],

and

r(x, y) = t(x,y) p(x, y)q(x,y).

If k = -1, let p(x,y) = x y zy2 + c2y, q(x,y), and r(x,y) as above.

Q.E.D.

Theorem 3.1 has the foolowing corollaries

COROLLARY 1 TO THEOREM 3.1. If

M n

t(x, y)= E E tijjxy

i=0 j=O

satisfies the symmetric or the skew symmetric property with respect to y, with too

and tmo nonzero, then there exist three polynomials p(x, y), q(x, y), and r(x, y) such

that

t(x,y) = p(x, y)q(x, y) + r(x, y),

where

2 2 m-2n-2

p(x,y) = pijxiyi, q(x, y)= Z qijxiyj,

i=o j=o i=0 j=o

and

m-1 n-1

r(x,y)= Y rijxiyj.

i=1 j=1

- 30

PROOF:

Note that since t(x, y) is either symmetric, or skew symmetric with respect to

y, ton, and tmn are also both nonzero. Following the notation used in the proof of

Theorem 3.1, the polynomials a(y), b(x), c(y), and d(x) are such that b(x) = d(x),

and a(y) and c(y) satisfy the symmetric property. Therefore, by Proposition 3.11

or Corollary 3.13, each root of a(y) and c(y) can be paired with its reciprocal, and

we can choose al(y),a2(y),cl(y), and c2(y) such that cala2 = 1, a3a4...On = 1,

7172 = 1, and 7374.---7n = 1. Thus, the conditions of Theorem 3.1 are satisfied, and

the result follows.

Q.E.D.

COROLLARY 2 TO THEOREM 3.1. If

n n

t(x, y)= Y E tijxyj

i=O j=o

satisfies the symmetric property with respect to both x and y and if the matrix

T = (tij) is symmetric, then there exist three polynomials p(x, y), q(x, y) and r(x, y)

such that :

t(x, y) = p(x, y)q(x, y) + r(x, y),

where

2 2 n-2n-2

Prx,y)= ^ pijpx'yy q(x,y)= qijx'y3

i=0 j=o i=0 j=O

and

n-ln-1

r(x,y) = rijx yj.

i=1 j=l

31 -

Moreover, the polynomials p(x,y),q(x,y), and r(x,y) can also be chosen to be sym-

metric with respect to both x and y.

PROOF:

In the case where too 0 0, the polynomials a(y), b(z), c(y) and d(x) defined

in Theorem 3.1 all have the same coefficients and therefore have the same decom-

positions. In the case, where too = ton = tm0 = tmnn = 0, the result follows from

Proposition 3.17.

Q.E.D.

3.4. The LU Factorization

For a detailed discussion of how the LU Factorization Theorem can be used to

implement operators in such a way that their operation count is reduced, see Nikias

et. al [21].

PROPOSITION 3.18. If a 5x5 matrix M has an LU factorization (without permu-

tations), then there exist two rank one matrices A and B and a matrix R = (rij)

such that

M= A+B+R,

and rlj = r2j = ril = ri2 = 0, for i,j = 1, ...5.

PROOF:

Suppose M can be written M = L U, where L is an lower triangular matrix,

and U is an upper triangular matrix. Suppose L = (lij), and U = (uij).

Let Ui be the 5x5 zero matrix whose ith row is the ith row of U and let

V = U3 + U4 + U5. M can then be rewritten as

32 -

M = L.U = L (Ui +U2 + V) = L- U1 + L -U2+L. V = A+B+R,

where

A = L U, B= L U2, and R= L V.

The matrix A is a rank one matrix of the form

A =

121

131

141

.151_

and the matrix B is a rank one matrix of the form

122

132

/42

.52.

[0 1 u23 u24 U25 .

On the other hand, R is a matrix whose first two rows and columns are all zero.

Q.E.D.

Thus, if t is a 5 x 5 operator, it can be viewed as a 5 x 5 matrix. If t has an

LU factorization without permutations, then there exist two rank one operators u,

and v and an operator r such that t = u + v + r.

Since rank one operators are separable (i.e. have an exact factorization), there

exist four 3 x 3 operators t1, t2, t3, t4, such that

u tlEDt2 and v=t3 f t4.

[I U12 U13 T114 U151

33 -

Therefore, t has the following decomposition

t = (tl t2) + (t3 t4) + r,

where

< tll > t12 t13

r = t21 t22 t23

t31 t32 t33 -

While Proposition 3.18 appears to be as good as the results in Corollary 3.2,

note that the center pixel in the image is at the (1,1) location. Thus, the data must

be shifted before it can be convolved with the operator r.

While the methods discussed in Sections 1-4 give techniques for the decom-

position of an n x n operator, most operators used in image processing are either

symmetric or skew symmetric with respect to either the x or the y axis. If n is an

odd integer, then the rank of the operator is no more than (n + 1)/2 in the symmet-

ric case, and no more than (n 1)/2 in the skew symmetric case. Using the rank

method or the LU factorization on symmetric or skew symmetric operators of low

rank, it is frequently the case that a decomposition can be produced which will be

nearly as efficient an implementation as that provided by our methods. Thus, the

methods presented here, are more apt to provide a decomposition superior to those

provided by the LU factorization on templates of relatively full rank.

3.5. Grobner Bases

Grobner bases have been introduced to Bruno Buchberger [4] in his disserta-

tion in 1965. They are named after W. Grobner [9], his thesis advisor. Buchberger's

34 -

method of Gr6bner bases is a technique that provides algorithmic solutions to a va-

riety of problems in ideal theory as well as solutions to numerous problems in many

other fields. It is a general purpose method for multivariate polynomial computa-

tions. It is implemented in all major computer algebra systems such as MACSYMA,

MAPLE, MuMATH, REDUCE, SCRATCHPAD II, etc... In spite of its simplicity,

the Buchberger algorithm solves a wide range of problems from symbolic algebra to

computational geometry.

In this dissertation, we will apply it to systems of algebraic equations and linear

homogeneous equations with polynomial coefficients.

3.5.1. Terminology

Let F be a field, and let F[yI, y2, ..., Yn] be the ring of n-variate polynomials over

F. The letters f, g, h will stand as typed variables for polynomials in F[yl, Y2,..., yn];

H, G for finite subsets of F[yl, Y2,..., yn]; t, u for power products of the form y'1 ...yi";

a,b for elements in the field F, and i,j for natural numbers.

DEFINITION 3.19. If H is defined as above, then the ideal generated by H is the set

Ideal(H)= {3 fhi : hi e H, fi E F[y, y2,..., ]}

of all linear combinations of the elements of H in F[y, y2, ..., Yn].

DEFINITION 3.20. A total ordering < on the power products y1 ...yn is called ad-

missible if 1 = y ..".y is minimal under < and multiplication by a power product

preserves the ordering.

35 -

Examples for such orderings are the total degree ordering (i.e. 1 < YI < Y2 <

y2 < YlY < Y2 < y3 < ... in the bivariate case), and the purely lexicographic

ordering (i.e. 1 < yl < 2 < ... < Y2 < y1Y2 < y2Y2 < ... < Y2 < YIY < ... in the

bivariate case).

Given a fixed ordering <, we will denote the coefficient of t in g by coeff(g, t), the

leading power product of g with respect to < by lpp(g), and the leading coefficient

of g with respect to < by Ic(g).

DEFINITION 3.21. Given a nonzero polynomial f, and two other polynomials g, and

h in F[yi,y2,...,yn], we say that g reduces to h modulo f or g --h f, if and only

if there exist u, b, b : 0, such that b.u.lpp(h) is identical with some monomial of g,

and f = g b.u.h.

Such a reduction can be viewed as a generalized division, deleting a monomial

in g. The result f is always strictly smaller than g, with respect to the ordering < .

Example 1. Let F = Q, g = y1y2y3 y2, and h = Y1y3 Y2. If we let u = y2, and

b = 1, and we use the lexicographic ordering, then f = g b.u.h = y2 y .

DEFINITION 3.22. We say that f is in normal form (or reduced form) if and only

if there is no f' such that f --h f.

Example 2. Let g and h be as in Example 1. If h1 = h, and h2 = y1Y2 Y1, then

f = y2 y2 of Example 1 is irreducible modulo {hl, h2}.

Another reduction of g is the following one: g -> f' = Y1Y3-Y2 = 2

In this case, h" is irreducible module {h, h

In this case, h" is irreducible module {hi, h2}.

36 -

Therefore, both f and f" are normal forms of g modulo {hl, h2}.

Hence, there exist several reduction paths that can lead to a normal form,

modulo a basis (or a set) H.

DEFINITION 3.23. The set H is called a Grobner basis if and only if each g has a

unique normal form modulo H.

We will note that different orderings will give rise to different Grobner bases,

and that a Gr6bner basis is not necessarily unique. It is possible to reduce a poly-

nomial of the basis modulo another polynomial of the basis, and hence obtain a

smaller basis.

As a matter of experience, using total degree ordering yields much better com-

puting time in general. On the other hand, when using the purely lexicographic

ordering, the resulting basis is in triangular form (i.e. each polynomial introduces

at most one new variable).

3.5.2. Application to Systems of Equations and Algorithm

Given a system of algebraic equations H, the Grobner basis can be used to:

-decide whether H is solvable;

-decide whether H has finitely or infinitely many solutions;

-find all solutions of H, in the finite case.

The following algorithm shows how the Gr6bner basis method using the lexico-

graphic ordering can be applied to determine the answer to a system of equations.

- 37 -

Algorithm

G = Grobner basis of H.

If 1 eG

then H is unsolvable.

Else if there exist an integer j such that no leading power product of the

polynomials in G is of the form yj, for some e,

then H has infinitely many solutions,

else do f= the (only) polynomial in G n F[yl],

X1 = {(a): f(a) = 0},

for i = 1,..., n 1,

Xi+1 = 0,

for all (al,..., ai) E Xi, do

K = {g(al,...,ai, Yi+l) :g G G n F[yl,..., yi+]

\F[yl,..., yi]},

f = g.c.d. of the polynomials in K,

Xi+i = Xi U {(a,...ai, a): f(a) = 0}.

Finally, Xn contains all the solutions of H.

The Gr6bner basis method was used in some examples of Chapter 5 to show

that a given template cannot be decomposed into a certain form. Several examples

are presented in Chapter 5. The equations were set up as a system which contains

both linear and nonlinear equations. For example, for a 5 x 5 template, we will

have 25 equations, 17 of which are not linear. To avoid having to use Gr6bner bases

theory to solve a system of nonlinear equations, we reduced the problem to roots

of polynomials of one variable to decompose the boundary of the template. Once

that outer boundary has been matched we only face a system of linear equations to

38 -

match the second boundary. While the Gr6bner bases method can be very slow in

solving nonlinear equations, it does a very fine job on solving linear equations. All

computations were done on MAPLE.

CHAPTER 4

VARIANT OPERATORS

We now turn our attention from results concerning shift-invariant templates

to analogous results for templates, which are not necessarily shift-invariant. Recall

that the Discrete Fourier Transform and the Gabor Transform are examples of vari-

ant templates. We will try to determine the class of templates that are separable

(i.e. decomposable into the convolution of a horizontal template and a vertical tem-

plate or vice-versa). Since not all templates are separable, we will show that every

template can be written as the sum of separable templates.

4.1. Decomposition of General Operators

DEFINITION 3.1. A template t defined on Z x Z has column rank one if it has finite

support and there exists a column vector u such that for each pixel location x in

Z x Z and every column vector v in t(x), there exists a scalar A such that v = Au.

A similar definition could be given for row rank one templates.

DEFINITION 4.2. A template t defined on Z x Z and having finite support is said

to be horizontal if, for every x in Z x Z, there exists a unique nonzero row vector v

39 -

40 -

in t(x). Similarly, a template t is said to be vertical if, for every x in Z x Z, there

exists a unique nonzero column vector v in t(x).

THEOREM 4.3. Ift is a column rank one template defined on Z x Z and has finite

support, then there exist an invariant vertical template t1 and a (possibly variant)

horizontal template t2 such that t = t1 6 t2.

PROOF:

For every x in Z x Z, t(x) is an m x n array. This array contains n column

vectors v1,V2,...,vn. Since the template has rank one, for each i = 1,...,n, there

exists a scalar Ai and a fixed column vector u such that vi = Aiu. Therefore,

t(x) = [u] [A1 A2 ... An] (outer product),

where [u] is an invariant vertical template and [A1 A2 ... An] is a (possibly variant)

horizontal template.

Q.E.D.

Remark 4.4. Let X be an n x n subset of Z x Z, and let F be a value set. The

two-dimensional Fourier template K is defined by

1

K(u,v) = -exp{-(uXi + vY1)27i/n},

n

where X1, Y1 C FX are defined by

Xi = {((x,y),z) : x = z} and YI = {((x,y),z) y =z}.

Even though it can be shown that K is a separable variant template, (i.e. it

splits into the product of a horizontal and a vertical), K has neither row rank one

-41 -

nor column rank one. If the Fourier template K did have column rank one, then by

Theorem 4.3, it could be factored into the product of an invariant vertical template

and a variant horizontal template.

DEFINITION 4.5. If a template t defined on Z x Z has finite support, then we say

that t has column rank at most k if there exist k column vectors ul,u2,..., uk such

that for every pixel location x in Z x Z and every column vector v in t(x), there

exist k scalars AI, A2,..., Ak such that v = k=1 Aiui.

Q.E.D.

THEOREM 4.6. If a column rank k template t defined on Z x Z has finite support,

then there exist k column rank one (possibly variant) templates tl,t2,...,tk such

that t = ti + t2 +... tk.

PROOF:

By Definition 4.5, for every pixel location x in Z x Z, and every column vector

vi in t(x), there exist k scalars Ail, Ai2,..., Aik such that vi= =1 Aijuj.

Let vij = Aijuj for i = 1,...,n. If tj(x) = [vj v2j ... vnj] for j = 1,...,k,

then each tj is a column rank one template and t = ti + t2 + ... + tk.

COROLLARY 4.7. If a template t defined on Z x Z has finite support and column

rank k, then there exist k vertical invariant templates ri, r2, ..., rk, and k horizontal

variant templates Sl,S2,..., sk, such that

k

t = re si.

i=1

- 42 -

PROOF:

This result follows immediately from Theorems 4.3 and 4.6.

Q.E.D.

Note that a result similar to Corollary 4.7 can be proved for row rank k tem-

plates.

THEOREM 4.8. Let t be a template on a finite rectangular subset of the plane with

m rows and n columns, and let Mt be the matrix described in section 4 of Chapter

2. If Ars denotes the matrix [mn(r1_)+i,n(j_-)+s] for i = 1,...,n and j = 1,...,n,

then the template t is separable if and only if each matrix Ars has rank one for all

r = 1,...,m and s = 1,...,m.

PROOF:

Since for every r = 1,...,m and s = 1,..., m, Ars is an n x n rank one matrix,

there exist a column vector Vrs and a row vector hrs such that Mrs = Vrshrs.

Let the column vectors V1s, 2s, ..., Vns form the sth n x n block of the block

diagonal matrix V, for s = 1, ..., m.

On the other hand, let the ith elements of the n row vectors hri, hr2, ..., hrn

form the diagonal elements of the (r, i) diagonal block of the block matrix H, for

r = 1, ...,m.

A direct computation V H will give the matrix Mt. Under the isomorphism b,

V and H represent the image of a vertical and a horizontal template, respectively.

We then conclude by Theorem 2.6, that t is separable.

Q.E.D.

43 -

COROLLARY 4.9. Let t be a template on a finite rectangular subset of the plane with

m rows and n columns, and let Mt be the matrix described in the previous paragraph.

If Ars denotes the matrix [mn(r_-)+i,n(j_-)+s] for i = 1,...,n and j = 1,...,n, then

the template t can be written as the sum of k separable templates if and only if each

matrix Ars has rank at most k for all r = 1,...,m and s = 1,...,m.

PROOF:

Since for every r = 1,...,m and s = 1,...,m, Ars is an n x n matrix of

at most rank k, there exist k column vectors v s,v2s,..., Vs and k row vectors

h1, h2.. hs. such that if As = vshs, for j 1,... k, then Ars = Al + A2s

... + As. For j = 1,..., k, Ars is an n x n rank one matrix, and therefore Theorem

4.9 applies. If for every j = 1,...,k Mj is the block matrix [A s], for r = 1,...,m

and s = 1, ..., m, then there exists two matrices Vj and Hj such that Ajs = Vj Hj.

Since Mt = Mll+M2 ...+Mk implies that Mt = VI- Hi+V2 H2-+...+VkHk,

Theorem 4.8 can now be used to conclude that the template t is the sum of k

separable templates.

Q.E.D.

Remark 4.10. Similarly, if we let Ars denote the matrix [mn(i-_)+r,n(s-l)+j], then

we can expect the same type of results as in Theorem 4.8. The only difference is that

the horizontal and the vertical templates are convolved in the reverse order. Recall

also that the two-dimensional Fourier template is separable. It is straightforward to

check that the Fourier template has the property that each of the matrices Ars has

rank one.

Remark 4.11. Note that according to the definition of the convolution of two tem-

44 -

plates t E s, it is as if we are "sliding" the first template t along the second one

s. Therefore, in order to minimize the number of equations and therefore minimize

the complexity of the solution, it will be very "helpful" to have the first template

t invariant. The goal of the next five theorems is to make this last statement more

precise. In fact, in this case, polynomials (rather than template notation) can be

used in the decomposition.

THEOREM 4.12. If t is a 5 x 5 variant template, then there exist seven templates

tl,t2,t3, t4,t5,t6, and t7 such that

t = (tl t2) + (t3 E t4) + (t5 ) t6) + t7,

where tl,t3, and t5 are 3x3 invariant templates and t2, t4,t6, and t7 are 3 x 3

variant templates.

PROOF:

Let t be of the form

t-2,-2 t-2,-1 t-2,0 t-2,1 t-2,2

t-1,-2 t-1,-1 t-1,0 t-1,1 t-1,2

t = t,-2 to,-1 < tO,0 > t0,1 to,2

tl,-2 tl,-1 t1,0 t,1 t1,2

t2,-2 t2,-1 t2,0 t2,1 t2,2

and let us define the following shift invariant templates:

- 45 -

0

<0>

0

0- 0

0 t3= 0

1 1

<0> 0 and t5=

0 0. .0

1

<0>

1

In addition, let us define the following variant templates:

t-2,-2

t-1,-2

0

t-2,-1

<0>

t2,1

t-2,0

<0>

t2,0

0

t1,2

t2,2-

t-2,1

<0>

t2,-1

t-2,2

t-1,2

0

0

t0,2

0

Note that a routine calculation using the definition of the operation can be

used to show that the template t7 = t (tl t2 + t3 E t4 + t5 D t6) is a 3 x 3

template. Thus, t can be decomposed as the sum and product of at most seven

3 x 3 templates.

Q.E.D.

tl=

and

0

t 1,-2

t2,-2

t6 t It,-2

46 -

4.2. Decomposition of Symmetric Operators

THEOREM 4.13. If t is a (2m + 1) x (2n + 1) variant template symmetric with

respect to both axes, then there exist four templates t1,t2,t3, and t4, such that

t = (tl t2) t3 + t4,

where t1 is a 3 x 3 invariant template, t2 is an (2m- 1) x (2n- 1) variant templates

symmetric with respect to both axes, t3 is a (2m + 1) x 1 vertical template, and t4

is a 1 x (2n + 1) horizontal template.

PROOF:

Let t2 be a (2m 1) x (2n 1) template symmetric with respect to both axes,

whose left top m x n corner is given by the equations:

St-m+i,-n+j if (i,j) =(0,0), (0,1), (1,0), (1,1)

t2 r(i,j) otherwise,

where

t-m+i,-n+j + t-m+i,-n+j-2 if i <

t-m+i,-n+j + t-m+i-2,-n+j if > j

r(i,j) =

t-m+i,-n+j + t-m+i-2,-n+j

+t-mn+i,-n+j-2 + t-m+i-2,-n+j-2 if i = j

Because of the symmetry of t2, it is easy to know all the other entries of that

template.

47-

1

Now if tl is the template t1 = 0

-1

0 -1-

< 0 > 0 then t (ti ( t2) is

0 1

simply the sum of a vertical template t3 and a horizontal template t4.

Q.E.D.

THEOREM 4.14. If t is a (2m + 1) x (2n + 1) variant template symmetric with

respect to both axes, then there exist five templates t1,t2,t3,t4, and t5 such that

t = (tl @ t2) + (t3 t4)+t5,

where t1 and t3 are 3 x 3 invariant templates and t2,t4, and t5 are (2m 1) x

(2n 1) variant templates symmetric with respect to both azes.

PROOF:

The main idea in this proof is to decompose the cruciform template s (i.e.

sij = 0, whenever i = 0 or j $ 0) derived from the proof of Theorem 4.14.

If

t 0 1 0>

t3= 1 < 0 > 1 ,

0 1 0

and t4 is the cruciform template of size (2m 1) x (2n 1) with same extremal

elements as in s, then a simple computation shows that t (ti E t2 + t3 E t4) is a

48 -

(2m + 1) x (2n + 1) template symmetric with respect to both axes.

Q.E.D.

THEOREM 4.15. If t is a (2m+ 1) x (21+ 1) variant template skew-symmetric with

respect to both the vertical and the horizontal axes, then there exist two templates tl

and t2 such that

t = tl t2,

where tl is the invariant template defined by

1 0 -1

t= 0 <0 > 0

-1 0 1

and t2 is a (2m 1) x (2n 1) variant template symmetric with respect to both

axes.

PROOF:

If t2 = [t2ij] is the (2m 1) x (2n 1) template symmetric with respect to

both axes defined in the proof of Theorem 4.13, then it is a matter of computation

to check that tl ( t2 is the given (2m + 1) x (2n + 1) skew-symmetric template t.

Q.E.D.

THEOREM 4.16. If t is a (2m + 1) x (2m + 1) totally symmetric template and if

the corner values are zero, then there exist three templates t1,t2, and t3 such that

t = (tl E t2) + t3,

49 -

where ti is a 3 x 3 invariant template and t2 and t3 are (2m 1) x (2m 1) totally

symmetric variant templates.

PROOF:

If the corner values of t are zero, then each boundary has either (2m 1) or

(2n 1) elements which form the boundaries of a (2m 1) x (2n 1) template t2

symmetric with respect to both axes.

0 1 0-

If ti = 1 < 0 > 1 then the template t3 = t ti 1 t2 is a (2m 1) x

.0 1 0.

(2n 1) template symmetric with respect to both axes.

Q.E.D.

CHAPTER 5

EXAMPLES

In the examples presented below, we have indicated the value in the center

pixel location by < x > Except for Example 6, all templates are shift-invariant.

Example 1.

The template

-13

2

7

2

-13

7

22

<27>

22

7

-13

2

7

2

-13

(used by Haralick [4]) satifies the conditions of Corollary 2 to Theorem 3.1. Using

the techniques of Corollary 2 it has the following decomposition

-19.73

< -28.9433 >

-19.73

-13 1

-13 1

-19.73 -.672

-13 1

50 -

-1.672

< 1.5411 >

-1.672

-13

-19.73

-13

1

-1.672 +[8.39]

1

- 51 -

Using the rank method (note that h has rank 2), the template h can be de-

composed as

S 1

h = 1.5282

1

+

1

1.2396

1

1.5275

< 2.3191 >

1.5275

1.3333

< 1.6528 >

1.3333

1 --13

1.5282 E 21.7364

1 -13

19.8575

< -33.2024 >

19.8575

1 2

1.2396 e < 14.5208 >

1 2

Using the LU factorization, the template h can be decomposed into

S1

h= 1.5182

1

1.5182

< 2.3049 >

1.5182

1 -13

1.5182 E 21.7364

1 -13

21.7364

< -36.3433 >

21.7364

[17.3077

+ 23.0769

.17.3077

23.0769

< 30.7692 >

23.0769

Thus, the first decomposition provides a more efficient implementation of h

than the other two decompositions.

-13

21.7364

-13

-13

21.7364

-13

17.3077

23.0769

17.3077

- 52 -

Example 2.

Using Corollary 3.2, the template

-5

-8

-10

-8

-5

-4

-10

-20

-10

-4

0

0

<0>

0

0

can be decomposed into the form

-4

< -10 >

-4

-5' 1 0 -1"

-8 D O <0> 0 + [-12 <0> 12].

-5 1 0 -1

Note that the template s is skew symmetric with respect to the y-axis and

symmetric with respect to the x-axis.

-5

s = -8

-5

- 53 -

Example 3.

The Laplace template

-1 0

0 0

<4> 0

0 0

-1 0

satifies the conditions of Proposition 3.6. It can be decomposed into

1 0 0- 0 0 -1-

1= 0 < 0 > 0 0 <0> 0 +[4].

0 0 1 -1 0 0

- 54 -

Example 4.

The template

01 0 10

00 0 00

t= 0 0 <0> 0 0

1 0 0 01

00 0 00

cannot be decomposed into the form (tlit2) + t3, where each ti is a 3 x 3 operator.

Using the LU factorization or the rank method, t can still be decomposed as the

sum and product of four 3 x 3 operators. Note that this template is symmetric with

respect to the y-axis. Thus, the hypothesis in Corollary 1 to Theorem 3.1 that the

corner values are different from zero is necessary.

- 55

Example 5.

The template

0

2

<1>

0

0

cannot be decomposed into the form (ti E t2) + (t3 t4), where each ti is a shift-

invariant 3 x 3 operator. Note that this template is symmetric with respect to the

y-axis.

- 56 -

Example 6.

Let X = {(1,1), (1,2), (2,1), (2,2)}, and let t be the variant template defined

by the rules

1 00 o0

t(l, 1) = t(1,2) = ,

S0 1 0

[1 0" 0 0"

t(2,1) = and t(2,2) -=

00 11

The matrix Mt associated with t is given by

Mt =

Note that while the template t has rank one for all x in X, it fails to be

separable.

CHAPTER 6

OPERATOR INVERSION

6.1. Introduction

Inverse problems have come to play a central role in modern applied mathemat-

ics, such as in mathematical physics, in imaging areas such as tomography, remote

sensing, and restoration. Rosenfeld and Kak [30] explained the equivalence between

techniques such as the Wiener filter and the least square methods and the problem

of inverting a block circulant matrix with circulant blocks. In 1973, G.E.Trapp [35]

showed how the Discrete Fourier Transform could be used to diagonalize and invert

a matrix that was either circulant or block circulant with circulant blocks. While the

algebraic relationship between circulant matrices and polynomials was completely

formulated by J.P. Davis [6], it was P.D.Gader and G.X.Ritter [9,27] who made

the connection between polynomials and circulant templates. A circulant template

defined on a rectangular array with m rows and n columns is invertible if and only

if its corresponding polynomial p(x,y) has the property that p(wA,w ) / 0, for

all 0 < j < n, and 0 < k < m. (The symbol wn denotes the root of the unity,

exp(27zi/n).)

An application of this method is that the usual 3 x 3 mean filter defined on a

rectangular array with m rows and n columns is invertible if and only if the number

-57 -

58 -

3 does not divide either m or n. Gader asked if a similar result applies to the von

Neumann mean filter. The von Neumann mean filter is the template t defined by

The main result of this chapter is to prove that if t is defined on a square array

with m rows and m columns, then t is invertible if and only if neither 5 nor 6 divide

m.

In this section, we will briefly describe the relationship between polynomial

algebra and image algebra. More precisely, the relationship between circulant tem-

plates and quotient rings of polynomial rings. Important results, such as how to

decompose circulant templates and how to obtain minimal local decompositions of

separable circulant templates can be found in Gader [8].

6.2. The von Neumann Template

6.2.1. Preliminaries

DEFINITION 6.1. A m x m matrix C = c(i,j) is said to be circulant if and only if

for every s E Z, cij = c(i+s)(mod m),(j+s)(mod m)-

Thus, C is a matrix of the form

Each row is just the previous row cycled forward

in each row are just a cyclic permutation of those in

C = circ(c, cl,...,c,_l1).

The m x m permutation matrix P =

one step, so that the entries

the first row. We will write

0

... 0

... 0

... 1

... 0

0 0 0

1 0 0

100

is called the basic

permutation matrix.

A circulant matrix C = circ(c, cl, ..., cm_1), can be associated with the poly-

nomial pc(x) = co + clx + ... + cm-1xm-1. If C = circ(co, cl, ..., cm-1), then

C = pc(P). Because of this representation, circulant matrices have a very nice

structure which can be related to P. Since a matrix is circulant if and only if

it commutes with the matrix P, the product of two circulants is again a circulant.

Thus, the set Cm of all m x m circulant matrices is a subring of the set of all matrices

of order m. Furthermore, circulant matrices commute under multiplication.

Let R[x] denote the ring of polynomials in one variable with coefficients in R,

and R[x]/(zm 1) the quotient ring of polynomials modulo zm 1. Multiplication

of elements of this ring are performed by replacing za" by 1.

- 59 -

co

Cm-1

c2

Cl

... Cm--

... Cm-2

60 -

If we define the mapping r7 : Cm -+ R[x]/(xm 1), by 77(C) = Pc(x), then 77

is a ring isomorphism [6].

DEFINITION 6.2. Let B be an mn x mn matrix. We say that B is block circulant

with circulant blocks of type (m, n) if and only if there exist n x n circulant matrices

CO, C1, ..., Cm-1 such that B = circ(Co, C1,..., Cm-1)

As in the circulant case, given a block circulant matrix of the form

C = circ(co, c, ..., cmn-1), with circulant blocks of type (m, n), we can associate

a polynomial pC(x, y) in two variables as follows:

pc(x,y) = (co + cly + ... + cnyn-l) + x(cn + Cn+1Y + ... + c2n-In-1) +

S+ xn-l(c(m-l)n + C(m-l)n+ly + ... Cmn-lyn-1).

As in the one variable case, R[x,y]/(xml 1,yn 1) denotes the quotient ring

of polynomials modulo xm 1, and yn 1.

If we define the mapping 0 : Cmn R[x,y]/(xm-, yn-1) by 0(C) = PC(, y),

then 0 is a ring isomorphism [6].

Therefore, given a circulant matrix M e Cmn, M is invertible if and only if

0(M) = pM(x,y) is a unit in R[x,y]/(xm 1,yn 1).

Let wn = exp(-27rij/m), where i = ---, and let A be the diagonal matrix

diag(wJ), for j = 0,1, ...,m 1. It can be proved that the permutation matrix

described above, satisfies the equation P = FAF*, where F is the one dimensional

Discrete Fourier Transform of order m, and denotes the conjugate transpose. Since

F is a unitary matrix, we have that F* = F-1. Therefore, if C is circulant, then

C = pc(P) = pc(FAF*) = Fpc(A)F* = FAF*,

where A is the diagonal matrix, diag(pc(wu)), for j = 0, 1, ..., m-1, whose diagonal

61 -

terms are the Fourier coefficients. Thus, any circulant matrix can be diagonalized

by the Fourier matrix.

Similarly, if F is the mn x mn two-dimensional Discrete Fourier matrix, and

M is a block circulant matrix with circulant blocks, then we can write M = FDF*,

where F is the mn x mn two-dimensional Fourier matrix, and D is a diagonal matrix.

Therefore, if the Fourier coefficients of M are all nonzero, then D is invertible

and the inverse of M is given by M-1 = FD-1F* [6].

Let t be an invariant circulant template. The following theorem has been

proved in Gader [8] .

THEOREM 6.3. Let t be a circulant template on the coordinate set X. The corre-

sponding matrix Mt is a block circulant matrix with circulant blocks.

Since we can write Mt = FDF*, where D = diag(pt(w, ,wn)) for

j = 0,...,m- 1, and k = 0, ..., n- 1, Mt is invertible if and only if pt(wn, wk) = 0,

for every j and k.

Example: Let t be the template defined on the m x n array X by

t(i,j) ={(i,j), (i + l(mod m),j), (i 1(mod m),j), (i,j + l(mod n))

(i,j 1(mod n))}.

t is called the von Neumann averaging (or mean) template. This template can also

be defined by the rule

S1 if i- k + j l < 1

S0 otherwise

62 -

The polynomial associated with this template, via the isomorphism 0, is

pt(x, y)= 1 + x + x"-1 + y + ym-1

Let X be a rectangular array with m rows and n columns, and let E be the

template defined by :

E(x) = {(y, ex(y)) : ex(x) = 1 and ex(y) = 0 if y = x}.

E is the identity template for the convolution D; that is

t E = E t= t.

A template t on X is said to be invertible with respect to E if there exists a

template s on X such that

t s=s t=E.

PROPOSITION 6.4. The von Neumann averaging template t, defined on a finite

rectangular coordinate set with m rows and n columns, is invertible if and only if

pt(wnm,w ) 0 for every j and k.

6.2.2. The Main Theorem

THEOREM 6.5. If X is an m x m rectangular coordinate set in Z x Z, and if t

denotes the von Neumann mean filter template defined on X, then t is invertible if

and only if neither 5 nor 6 divide m.

- 63 -

PROOF:

For j, k = 0,..., m 1 we have

Pt(wiw ) = 1 + w? + w- + wk + k

= 1 + 2cos(27rj/m) + 2cos(27k/m).

For convenience, we will write w instead of wi.

We first assume that one of the integers 5 or 6 divides m.

Let m = 5d, for some d. If we let j = d and k = 2d, then

pt(wJ,wk) = 1 + wd + -d 2d -+ -2d

-1 2 -2

= 1 + w5 + +5 + W5

2 3 4

= 1 +w swg + + 5 = 0

for this is the sum of the five fifth roots of unity.

Let m = 6d, for some d. If we let j = d and k = 3d, then

pt(wk) = +wd + -d +w3d w-3d

= 1 + 2cos(7) + 2cos(7/3) = 0.

Hence if 5 or 6 divides m, there exist two integers j and k such that pt( m, m)

= 0 and, by Proposition 2.4, the template t is not invertible.

WE now prove the sufficient condition.

Due to the symmetry of j and k, in the polynomial expression pt (w, w), it is

enough to restrict these integers to 0 < j < k < m/2.

Let j be the minimal value such that pt(wj,wk) = 0.

Note that if

q(x) = 1 + x + x(m-1)j + xk + X(m-1)k,

then the equation pt(wJ,wk) = 0, is equivalent to q(w) = 0. If we use the fact that

cyclotomic polynomials are the minimal polynomials of the primitive roots of unity

[3], then we have that q(ws) = 0, for any s relatively prime to m.

64 -

Claim: If pt(w,wk) = 0, then m/6 < j < k < m/2.

If for any z = r*exp(ia), Arg(z) = a, then Arg(wJ) = Arg(w-J). Note that

in order to find two integers j and k such that pt(w m,w') = 0, it is necessary

that r/3 < Arg(w-) _< 7. Since Arg(wj) = 2rj/mn, the previous argument gives

m/6 < j < m/2. The claim is proved.

We will consider two cases: j divides m and j does not divide m.

Case 1: j divides m.

Using the claim, we know that 2j < m < 6j. If j divides m, then m is of the

form 2j, 3j, 4j, 5j or 6j.

If m = 2j, then pt(wj,wk) = Oimplies that cos(2k/mn) = -1/2, which implies

that k = m/3. Since k is an integer, 3 divides m. Hence 6 divides m.

If m = 3j, then pt(wJ,wk) = 0 implies that cos(27k/m) = 0, which implies

that k = m/6. Since 3 and 4 divide m, 6 divides m.

If m = 4j, then pt(wj,wk) = 0 implies that cos(27rk/m) = -1/2, which implies

that k = m/3. Since 3 and 4 divide m, 6 divides m.

If m = 5j or 6j, then since j is an integer, we conclude that 5 or 6 divides m.

Case 2: j does not divide m.

Using the claim again, we can write m as 2j + r, 3j + r, 4j + r, or 5j + r for

some nonzero r < j.

If m = 3j + r, then either 2 divides m or 2 does not divide m. In the case

where 2 divides m, if 3 divides m, we are done. Otherwise, the greatest

commun divisor of m and 3, denoted by gcd(m,3), is equal to 1. Thus

q(w3) = 1 +w3j + j +W3k +3+ -3k = 0.

But 3j = m r, thus,

65 -

q(w3) = 1 + wm- + w-(m-r) + w3k + w-3k

= 1 + wr + w- + w3k + -3k,

which contradicts the minimality ofj. In the case where 2 does not divide

m, then gcd(m,4)=1, and

q(w4) = 1 +w4 + -4j + 4k + -4k = 0.

But, 4j = m r + j, thus,

q(w4) = 1 + m-r+j + -(m-r+j) + 4k + w-4k

= 1 + j-r -(j-r) 4k -4k,

which contradicts the minimality of j.

If m = 4j + r, then either 3 divides m or 3 does not divide m. In the case where

3 divides m, if 2 divides m we are done. Otherwise, gcd(m,4)=l. Thus,

q(w4) = 1 + w4j + w-4 + w4k + -4k = 0.

But, 4j = m r, thus,

q(w4) = 1 + wn-r + w-(m-) + w4k + -4k

1 +wr+ -r + w5k + -5k,

which contradicts the minimality of J. In the case where 3 does not divide

m, then gcd(m,9)=1 and

q(w9) = 1 + w9 + w-9j 9k + + -9k = 0.

But 9j = 2m 2r + j, thus,

q(w9) 1 + w2m7-2r+j + -(2m1-2r+j) + 9k + w -9k

= 1 + -2r + -(j-2r) + 9k + -,

which contradicts the minimality of j, since 1j 2rl < J.

If m = 5j + r, then if 5 divides m we are done. Otherwise, gcd(m,5) = 1 and

q(w5) = 1 +w5j + -5j +w5k + w-5 = 0,

since w5 is a primitive mth root of unity. But 5j = m r, thus,

66 -

q(w5) = 1wm + c-r -(m-r) + 5k + w-5k

1 + wr + w-r + w5k + -5k,

which contradicts the minimality of j, since r <_ j. This concludes

the proof of Theorem 3.1

Q.E.D.

A generalization of the von Neumann averaging template is the one defined by

t a if r = i, and s =

b otherwise

b

t = b b

b

where a, and b are two real numbers.

The associated polynomial is p(x, y) = a + bx + bxm-1 + by + bym-l, and

p(w ,wk) = a + bwj + bw-j + k + b-k

= a + b(wi + w-j +k + -k).

PROPOSITION 6.6. If a = 0, then the template t is invertible if and only if m is

an odd number.

PROOF:

Of course, we will assume, here, that b 7 0, in which case it can be factored.

If p(wj, k) = wj w+ w-j + wk + c-k

= 2cos(27rj/m) + 2cos(27rk/m) = 0,

-67 -

then the equation is equivalent to cos(27rj/m) = -cos(27rk/m), which implies

27rj/m + 27rk/m = 7r, and therefore, j + k = m/2. Since j and k are integers,

this proves that for any j, and k = m/2 -j, p(wj,wk) = 0. Therefore, the template

t is not invertible if m is an even number.

Q.E.D.

PROPOSITION 6.7. If a $ 0 and la/b| > 4, then the template t is invertible.

PROOF:

If p(wj,wk) = a + + b bw- + bw + bwk-1

= a + 2b(cos(2nj/m) + 2cos(27rk/m)) = 0,

then Icos(27rj/m) + cos(27rk/m)l = ja/2bI > 2, which is a contradiction. Therefore,

p(wJ,wk) is always nonzero, and the template t is invertible.

Q.E.D.

6.3. The Characteristic Polynomial Method

Let X be a coordinate set with m rows and n columns, and let t be a template

defined on X. Using the relationship between the matrix algebra and the image

algebra, under the isomorphism 0, we know that the corresponding matrix 0(t) =

Mt is an mn x mn block matrix. There are m2 blocks and each block is an n x n

matrix.

The determinant det(Mt AImn), where Imn is the unit mn x mn unit or

identity matrix, is called the characteristic polynomial of Mt. It has the form

p(A) = Amn + amn-_imn-1 + ... + alA + ao.

68 -

In p(A), the constant a0 is the determinant of Mt, while the coefficient amn-1

is the trace of Mt [12].

THEOREM 6.8. If ao is nonzero, then Mt is invertible and its inverse is given by

the formula:

M 1 = (Mmn-1 + a, 1Mn-2 + ... + allmn).

PROOF:

By the Cayley-Hamilton theorem, we know that the matrix Mt satisfies its

characteristic polynomial p(A), i.e p(Mt) = 0. Therefore,

M"nn + amn-iMmn-1 + ... + aiMt + aolmn = 0.

The previous equation can be rewritten as

Mtmn + amnlMtn-1 + ... + aMt = -aoImn, which implies

Mt[o (Mt- + amn-1tA -2 + .. alInn)] = Inn.

This last equation gives

Mt1 = M(ltmn- 1 + an-M~ n-2 + ... + al/mn).

Q.E.D.

Let E be the identity template defined in Section 6.2.1, and let tn be defined

as t t E ... t n times.

69 -

From the previous theorem, we can conclude that, given a shift invariant tem-

plate t, we can write its inverse as

-1 1

t_i =-(t1mn-1 + amn-2tmn-2 + ... + alE).

a0

Since the inverse of a block toeplitz matrix with toeplitz blocks is not neces-

sarily a matrix with the same properties, we see that even the inverse of a very

small size shift invariant operator can be (and most probably will be) a full size

variant operator. The previous theorem is one way to write the inverse of a shift

invariant template in terms of the template itself and therefore, in terms of invariant

templates.

6.4. Inverse of General Shift Invariant Operators

In this section, we discuss some properties of block Toeplitz matrices with

Toeplitz blocks as they are directly related to shift invariant templates which are,

as mentioned previously, used to implement convolutions. The idea here is to apply

some results on Toeplitz matrices to the image algebra.

As in the beginning of this chapter, the matrices of this section are indexed by

sets of the form {0, 1,..., m 1}.

DEFINITION 6.9. Let M = (mj) be an in x rn matrix. We say that M is a

Toeplitz matrix if and only if for every i,j E {0, 1,..., m 1}, and k in Z such that

i + k,j + k {0,1, ...,m 1}, we have aij = ai+k,j+k.

a0 a1 ... arm-

a-1 ao ... am_2

Thus, M is a matrix of the form

-a-(m-l) a-(m-2) ... a0

- 70 -

A Toeplitz matrix is constant along any diagonal and has its ijth entry as a

function of (i j) rather than of i and j separately.

DEFINITION 6.10. Let M = (Mij) be an m x m

an n x n matrix. We say that M is block Toeplitz

M is of the form

Ao

A_1

LA_(?m-) A-(m2)

where for every i in {-(m 1), (m 1)}, each Ai

block matrix, where each Mij is

with Toeplitz blocks if and only if

... Am-1

... Am_2

... Ao -

is a Toeplitz matrix.

The following theorem was proved in Gader [8].

THEOREM 6.11. Let X be a rectangular coordinate set. If t is a template on X,

and Mt is the corresponding matrix under 0b, then Mt is a block Toeplitz matrix

with Toeplitz blocks.

Since the inverse of a Toeplitz matrix is not necessarily a matrix of the same

form, we see that this very special structure of Toeplitz and block Toeplitz matrices

strongly suggests that an inversion scheme exploiting this structure would yield

significant savings in time and effort. Several researchers have tried to investigate

Toeplitz matrices and their inverses. Some algorithms like the Trench algorithm

proved to be powerful [33]. Other methods yielded very elegant formulas for the

- 71 -

inverse [1].

Let us first give some facts pertinent only to Toeplitz matrices.

Let M be an m x m Toeplitz matrix; u, v, two mxl vectors; and ei the transpose

of the vector (0, ...0, 1,0, ...0), where 1 is at the ith location.

Fact 1: If AM(o, ui, ..., um-l)t = (vo, Vl, ..., vm-1)t

then (um-1, Um-2, ..., uo0)M = (Vm-1, Vm-2, ..., vO)

Fact 2: If Mu = ep, and Mv = eq,

then Um-1-q = Vm-1-p.

Note that we have

[(Vm-l, ..., vo)M](uo, ..., tr-l)t = (m-,..., -, v)[M(uo,..., Um-1)]1

We will call a block vector of dimension m, a vector in which all entries are

n x n matrices. It will be denoted by a bold letter. We will denote by ei the vector

(0, ..., 0, I, 0, ..., 0), where the n x n unit or identity matrix I is at the ith entry, and

where 0 stands for the n x n zero matrix. L(Ao, A, ..., Am-1) represents the lower

triangular block Toeplitz matrix

Ao 0 ... 0

AI Ao ... 0

L=

-A(m-) A-(m-2) ... Ao

72 -

In a manner similar to the lower triangular matrix, we denote the upper trian-

gular block Toeplitz matrix by U(Ao, A, ..., Am-1) and finally, we let the matrix S

denote the shift matrix S = L(0, I, 0,..., 0).

PROPOSITION 6.12. Let M be an m x m block matrix with n x n blocks. If there

exist column block vectors x, and y such that

Mx = eo and My = SMe,,,

then M is invertible. Moreover, let w and z be two block vectors such that

wM = et and zMy = et)MS,

then

M-1 = A1B1 A2B2,

where

A1 = L(xo,xl,...,xnl-) ,A2 = L(-yo, -Yi, ...,-ym-),

B1 = U(I,-z0,...,-Zm_2) ,B2 = U(0, w,...,wm_2).

In the implementation of this theorem, we transpose the dreadful problem of

inverting an mn x mn block Toeplitz matrix with Toeplitz blocks to the problem of

solving 4n linear systems of mn equations with mn variables.

73 -

In the application of this theorem to image algebra, the result is quite inter-

esting. If t is an invariant template satisfying the conditions of the theorem, then

there exist four templates t1,t2,t3, and t4 such that

t-1 = tl t2 + t t4,

where t1, and t3 are shift invariant along each row, and t2 and t4 are shift invariant

along each column

CHAPTER 7

CONCLUSION AND SUGGESTIONS FOR FURTHER RESEARCH

This dissertation has addressed the problems of template decomposition and

template inversion. The framework is based on the image algebra and its relationship

to polynomial and matrix algebra.

In particular, it was shown that:

1) A two-dimensional shift invariant operators defined on a rectangular array

can be decomposed into sums and products of smaller size operators. The

decomposition methods presented here, are suitable for a machine with a

pipeline architecture.

2) Two-dimensional shift invariant operators exhibiting some kind of symme-

try require a much smaller number of operators in the decomposition. We

focused our attention on this type of operator because of their frequent use

in image processing.

3) The techniques were extended to those operators which may not be shift

invariant. Necessary and sufficient conditions were proved to test for the

separability of a variant operator. Necessary and sufficient conditions were

proved to test whether or not a variant operator can be written as the sum

of separable operators.

4) Necessary and sufficient conditions were given to characterize when the mean

filter based on the von Neumann configuration is invertible. Since the inverse

74 -

75 -

of a shift invariant operator is not necessarily shift invariant (actually, there

is very little chance it is shift invariant), methods were provided for writing

the inverse in terms of shift invariant operators.

We now list some suggestions for further research.

First, generalize the results obtained in Chapter 3 for two-dimensional shift

invariant operators to arbitrary shift invariant operators of different configuration.

Circular and convex configurations would be of particular interest. Can the poly-

nomial decomposition method be generalized to these more general operators ?

The decomposition methods for variant operators would certainly be imple-

mented if the cost in computation time was not prohibitive. Is there a better and

more direct approach to their decomposiiton ? Many methods given in Chapter

4 involve shift invariant and variant operators in the decomposition of variant op-

erators. A method involving only shift invariant operators will surely be a more

efficient one. An important example of a variant template is the Gabor Transform.

Can the results presented here be used to decompose this operator.

A generalization of the condition of the invertibility of the von Neumann mean

filter defined on a rectangular array, rather than on a square array should be proved

in the the very near future. The inverse of a shift invariant operator can be a full size

variant operator. Methods to write this inverse in terms of an acceptable number

of shift invariant operators can find their way in areas such as tomography, remote

sensing, and cryptography among others.

REFERENCES

1. A. Ben-Artzi and T. Shalom, "On the Inverse of Toeplitz and Close to Toeplitz

Matrices", Linear Algebra and Applications 75 (1986), pp. 171-192.

2. G. Birkhoff and J. Lipson, "Heterogeneous Algebras", J. Combinatorial Theory

8 (1970), pp. 115-133.

3. J. P. Blahut, "Fast Algorithms for Digital Signal Processing," John Wiley and

Sons, New York, 1985.

4. P. Buchberger, "Grobner Bases: An Algorithmic Method in Polynomial Theory",

Recent Trends in Multidimensional Systems Theory (1985).

5. J. L. Davidson, "Lattice Structure in the Image Algebra and their Applications"

to image processing", Ph.D. Dissertation, University of Florida, Gainesville, Fl

(1989).

6. J. P. Davis, "Circulant Matrices," John Wiley and Sons, New York, 1979.

7. D. Gabor, "Theory of Communications", J. Inst. Elec. Eng. (London) 93

(1946), pp. 149-157.

8. P. D. Gader, "Image Algebra Techniques for Parallel Computation of Discrete

Fourier Transforms and General Linear Transforms", Ph.D. Dissertation, Uni-

versity of Florida, Gainesville, Fl (1986).

9. W. Grobner, "Modern Algebraic Geometry," Springier Verlag, 1949.

10. R. M. Haralick, "Rlidges and Valleys on Digital Images", Journal of Computer

Vision, Graphics and Image Processing 22 (1983), pp. 28-38.

11. C. E. Heil and D. F. Walnut, "Continuous and Discrete Wavelets Transforms",

SIAM Review 31 No 4 (1989), pp. 628-666.

12. N. Jacobson, "Basic Algebra I," Freeman and Company, San Francisco, 1974.

13. D. E. Knuth, "The Art of Computer Programming, Sorting and Searching,"

Addison Wesley, Reading, MA, 1973.

14. D. E. Knuth and P. B. Bendix, "Simple Word Problems in Universal Algebra",

Proc. Computational Problems in Abstract Algebra, Oxford, England (1970),

pp.263-298.

- 76 -

- 77 -

15. B. Kutzler and S. Stifter (Grobner), "Automated Geometry Theorems Prov-

ing Using Buchberger's Algorithm", Proceedings SYMSAC '86, Wateloo,Canada

(1986, ACM publications), pp.209-214.

16. D. Li, "Recursive Operations in Image Algebra and their Applicaitons to Image

Processing", Ph.D. Dissertation, University of Florida, Gainesville, Fl (1990).

17. Z. M. Manseur and D. C. Wilson, "Decomposition Methods for Convolution

Operators", submitted.

18. G. Matheron, "Development of a Mathematical Structure for Image Processing",

Technical report, Optical Division Tech. Report, Perkin-Elmer (1983).

19. P. E. Miller, "Random Sets and Integral Geometry," Wiley, New York, 1975.

20. D. R. Musser, "Multivariate Polynomial Factorization", Journal of the ACM

22(2) (April, 1975), pp. 291-308.

21. C. L. Nikias, A. P. Chrysafis, and A. N. Venetsanopoulos, "The LU Decom-

position Theorem and its Implications to the Realization of Two-Dimensional

Digital Filters", IEEE Transactions on Acoustics, Speech and Signal Processing

3 (1985), pp. 694-711.

22. D. P. O'Leary, "Some Algorithms for Approximating Convolutions", Journal of

Computer Vision, Graphics and Image Processing 41 (1988), pp. 333-345.

23. G. X. Ritter, "An Algebra of Images", Proc. Conf on Intelligent Systems and

Machines (1984), pp. 333-339.

24. G. X. Ritter, "The Language of Massively Parallel Processing Computers", Proc.

1984 S. E. Reg. ACM Conference, Friendly Systems 1984-2000? (1984), pp.

224-231.

25. G. X. Ritter and P. D. Gader, "Image Algebra Implementation on Cellular Array

Computers", IEEE Computer Society Workshop on Computer Architecture for

Pattern Analysis and Image Database Management (1985), pp. 430-438.

26. G. X. Ritter, P. D. Gader, and J. L. Davidson, "Automated Bridge Detection

in FLIR Images ", In Proceedings of the Eighth International Conference on

Pattern Recognition. Paris, France (to appear).

27. G. X. Ritter and P. D. Gader, "Image Algebra Techniques for Parallel Image

Processing", Journal of Parallel Distrib. Comput. 4, No 5 (1987), pp. 7-44.

28. G. X. Ritter, M. A. Shrader-Frechette, and J.N. Wilson, "Image Algebra: A

Rigorous and Translucent Way of Expressing all Image Processing Operations",

In Proc. of the 1987 SPIE Tech. Symp. Southeast on Optics, Elec.-Opt., and

Sensors. Orlando, FL, (May 1987), pp. 7-44.

29. G. X. Ritter, J. N. Wilson, and J. L. Davidson, "Image algebra: An overview",

Journal of Computer Vision, Graphics and Image processing 49 (1990), pp.

297-331.

- 78 -

30. A. Rosenfeld and A. C. Kak, "Digital Picture Processing," Academic Press, New

York, London, 1982.

31. P. C. Sabatier, "Inverse Problems: An Interdisciplinary Study," Academic Press,

New York, 1987, pp. pp. 21-67, 99-175, 217-224.

32. J. Serra,, "Image Analysis," Academic Press, London, 1982.

33. S. R. Sternberg, "Language and Architecture for Parallel Image Processing",

Proc. of the conf. on Pattern Recognition in Practice, Amsterdam (1980).

34. S. R. Sternberg, "Biomedical Image Processing", Computer 16(1) (1983), pp.

22-34.

35. S. R. Sternberg, "Overview of Image Algebra and Related Issues", Integrated

Technology for Parallel Image Processing, Academic Press, London (1985).

36. G. E. Trapp, "Inverses of Circulant Matrices and Block Circulant Matrices",

Kyungpook Math. J. 13 No 1 (1973), pp. 11-20.

37. W. F. Trench, "An Algorithm for the Inversion of Finite Toeplitz Matrices",

Jour. SIAM 12,3 (1964), pp. 515-522.

38. J. N. Wilson, D. C. Wilson and G. X. Ritter, "The Image Algebra Fortran

Language", UF-CIS Technical Report TR-88-03,Center for Computer Vision

Research, University of Florida (Florida, June 1988).

BIOGRAPHICAL SKETCH

Zohra Z. Manseur was born in Ouenza, eastern Algeria, shortly after the be-

ginning of the Algerian war against France and grew up in Algiers. She received

a Licence-es-Sciences of Teaching in Mathematics from the University of Algiers in

1981, and a Diplome d'Etudes Superieures in Algebra in 1983 from the University of

Science and Technology of Algiers. She joined the graduate program in Mathemat-

ics at the University of Florida in August 1984 and received a Master's degree in

August 1986. She married Rachid Manseur in 1978 and has a son and a daughter.

- 79 -

I certify that I have read this study and that in my opinion it conforms to

acceptable standards of scholarly presentation and is fully adequate, in scope and

quality, as a dissertation for the degree of Doctor of Philosophy.

Ger rd X. Ritter, Chairman

professor of Mathematics

Professor of Computer and Information Sciences

I certify that I have read this study and that in my opinion it conforms to

acceptable standards of scholarly presentation and is fully adequate,in scope and

quality, as a dissertation for the degree of Doctor of Philosophy.

David C. Wilson, CoChairman

Professor of Mathematics

I certify that I have read this study and that in my opinion it conforms to

acceptable standards of scholarly presentation and is fully adequate, in scope and

quality, as a dissertation for the degree of Doctor of Philosophy.

ames E. Keesling

SPressor of Mathematics

I certify that I have read this study and that in my opinion it conforms to

acceptable standards of scholarly presentation and is fully adequate, in scope and

quality, as a dissertation for the degree of Doctor of Philosophy.

Li-Chien Shen

Associate Professor of Mathematics

I certify that I have read this study and that in my opinion it conforms to

acceptable standards of scholarly presentation and is fully adequate, in scope and

quality, as a dissertation for the degree of Doctor of Philosophy.

Soseh1 N. Wilson

Assistant Professor of Computer and Information

Sciences

This dissertation was submitted to the Graduate Faculty of the Department of

Mathematics in the College of Liberal Arts and Sciences and to Graduate School

and was accepted as partial fulfillment of the requirements for the degree of Doctor

of Philosophy.

December 1990

Dean, Graduate School