Matrix newton interpolation and progressive 3D imaging: PC-based computation

MATHEMATICAL COMPUTER MODELLING PERGAMON

Mathematical

and Computer

Modelling

35 (2002) 303-322 www.elsevier.com/locate/mcm

Matrix Newton Interpolation and Progressive 3D Imaging: PC-Based Computation E. DEFEZ Departamento

de MatemBtica

Aplicada,

Universidacl

46071 Valencia, A. Department

of Computer Waterloo,

J.

VILLANUEVA-OLLER

Departamento

de Matematica

Politecnica

de Valencia

Spain

LAW Scence,

Ontario,

University

Canada AND

Aplicada,

R.

J.

Universidad

46071 Valencia,

of Waterloo

N2L 3Gl VILLANUEVA’ Politecnica

de Valencia

Spain

rjvillanQmat.upv.es

(Received April 2001; accepted May 2001)

Abstract-For polynomials P(r) = A,?’ + An_~Pel + + Ala: + A0 in a real scalar z, but with coefficients Aj that are rectangular matrices, a generalization of Newton’s divided difference interpolatory scheme is developed. Instances of P(z) at nodes ri may be interpreted as slices of a digital 3D object. Mathematics code for this machinery is given and its effectiveness illustrated for progressively-transmitted renderings. Analysis, with supporting Mathematics code, is extended to a piecewise matrix polynomial situation, to produce practicable software for a PC-based computational system. Two experiments about 3D progressive imaging, employing a 6 Mbyte data base consisting of 93 CT slices of a human head, are discussed along with PC-based performance evaluation. How a 3D object is decomposed into 2D subsets in preparation for progressive transmission, as well as their selected ordering for transmission, are seen to affect quality of the emerging reconstructions, Extension to 4D objects is also discussed briefly, to provide introduction to, for example, application of matrix polynomial machinery within the field of functional magnetic resonance imaging. @ 2002 Elsevier Science Ltd. All rights reserved. Keywords-Progressive reconstruction, PC-based

transmission of images, Matrix Newton interpolation, progressive rendering.

Matrix polynomial

1. INTRODUCTION Modern medical imaging began with Roentgen’s discovery of X-rays, just over 100 years ago, and much of the rapid development in the science has occurred in the last 25 years. As the technology has matured, the quality and quantity of data produced in a typical examination have increased dramatically. For example, a computerized tomography (CT) exam of the early 1970s usually consisted of one or two slices and perhaps 100 kilobytes of data. By the early 1980s exams This research has been partially supported by the Plan de Incentive a la Investigation/2000 Politkcnica de Valencia. *Author to whom all correspondence should be addressed. 0895-7177/02/% - see front matter @ 2002 Elsevier Science Ltd. All rights reserved. PII: SO895-7177(01)00167-4

Typeset

de la Universidad

by &&-TI$

304

E. DEFEZ et al.

contained lo-20 slices and l-2 megabytes (MB) of data, and by the end of the decade, 50 or more slices with 25 MB of data per exam were common. Current spiral-CT scanners can image hundreds of slices per exam and produce 100 MB or more of data, while magnetic resonance or 3D ultrasound imaging is capable of capturing 509 MB of data in just a few minutes. Other imaging fields have also seen burgeoning amounts of data, as user requirements become more sophisticated.

In general, interactive video, for example, even

picture browsing drains computational and storage resources, while developments for full interactivity currently remain in the realm of virtual reality technologies research and development, The underlying major problem, a diverse area of active research, is how to manipulate these ever-increasing amounts of medical, or other, data.

Their sheer quantity creates problems in transmitting, viewing, interacting with, and storing results of processing. A common approach attempted is to compress the image data, somehow, and schemes under consideration typically try to create mathematical models of information in the image, and then store only “significant” portions of these models, thereby reducing the total amount of data. In such lossy schemes, the greater the compression, the greater the amount of information discarded, hence, generally, the poorer the quality of the compressed-then-decompressed image, for subsequent analysis. In medical, as well as other application areas of imaging science, differing approaches and techniques are under investigation for compression, transmission, and subsequent decompression of 3D images. One class of techniques imposes some compression scheme on the full digital object, hence, the transmitted image is not rendered until all transmission is finished. The traditional standard here has been a joint photographic experts group (JPEG) algorithm, based in part on the discrete cosine transform and designed for low compression rates. Recently, several algorithms for compression based, for example, on the wavelet transform have also been considered [1,2]. A second approach for developing techniques, and this is the scenario applicable for the development below, uses instead, progressive transmission for subsequent viewing of the 3D object. With such progressive rendering, a high degree of effective compression is achieved by not transmitting the remainder of the data in the event that the image is appearing to be not of interest. Or, during progressive viewing, if some feature is emerging as important, the receiver can interact with the transmission process and request contiguous subdata (but still lossless) immediately. Some mathematical machinery for progressive transmission has been investigated recently using decorrelation techniques [3], or wavelets [4], and more recently, emerging mathematical structure of matrix orthogonal polynomials [5], for example. Other approaches may be found in [6,7]. The first purpose of this paper is to introduce interpolation for matrix polynomials, and develop a Newton divided difference formulation that indeed generalizes such a structure from the familiar scalar polynomial situation. Then this matrix polynomial interpolatory machinery is used in a process for progressive transmission of a digital 3D object, composed of an ordered sequence of 2D subsets. Effectiveness of the theory is illustrated through application to a standard data set which has been employed recently for illustrating some other progressive imaging algorithms [5,8-lo]. The set in initial form represents 93 parallel CT slices of a human head, the slices each being of resolution 256 x 256. Not unsurprisingly perhaps, this straightforward implementation can consume considerable processing time for its evolving 3D renderings. To avoid the major computational overhead associated with this higher-degreed matrix polynomial scheme, the mathematical underpinning is extended to piecewise-polynomial structures, and linear piecewise segments are found to be (a) effective visually for progressive transmission/rendering, as well as (b) practicable for a single processor PC-based implementation, for 3D objects. Features of this extended scheme are illustrated through results from a first experiment. With the tool of an acceptable PC-based implementation (in Mathematics [ll]) of the piecewiselinear algorithm available, another major point, included in this paper, could be addressed empirically, and a second experiment exploits the same test CT data set, but reconfigured as a

Matrix Newton Interpolation

305

stack of 256 slices, each therefore containing 93 × 256 pixels. Both computational performance and effectiveness of the imaging are seen to be different between this situation and the first one. This addresses a basic point for progressive imaging: how the subsets of a 3D object are selected initially, as well as the ordering scheme imposed for carrying on with their transmission, can be crucial. The paper is organized as follows. Section 2 develops matrix Newton interpolation through generalization of divided differences from the scalar-polynomial situation. In Section 3, the matrix polynomial interpolant for progressive transmission of 3D images is discussed, with its visual effectiveness and computational costs illustrated for one section of the head d a t a set. In Section 4, notation for progressive transmission of an image is meshed with t h a t for piecewise (linear) matrix Newton interpolation, and error expressions are provided for numerically assessing goodness of corresponding reconstructions. Full Mathematica code for the computational schemes is given in Section 5. Two experiments are developed in Section 6, for progressive transmission/renderings generated from the same 3D image, but with their 2D subset sequences structured differently. Section 7 contains a general discussion of findings and emerging questions for future work. In the field of progressive transmission, a 3D digital object is usually stored as, in essence, a stack of parallel 2D slices of it--i.e., as an ordered collection of matrices with their common size corresponding to image resolution, say r × s. The ordering scalar could be, e.g., the distance along the axis perpendicular to the slices. Thus, progressive transmission of a 3D object NI = {(xi, Mi)}n=o, where x0 < xl < x2 < . " < Xn are distinct real numbers and the ~li are given r x s matrices (CT or MRI slices, for instance), can be accomplished with the following steps. 1 ~li 1,~lnl In the first step, the sender extracts a subset NIl C NI, NIl {( :;Ci, JJi=l, and transmits its elements. The recipient receives Nil and performs some reconstruction of NI. In the second step, the sender extracts a subset NI2 C NI ~ M1, to avoid sending the same information twice• 2 M i2~ln2 Therefore, the recipient has received NIl UNI2 = {( x i, JJi=l images, and again, he reconstructs an approximation to NI. In step k, the sender extracts a subset NIk C NI ~ (NIl U ... U N[k-1) to 1/-k _~lk~u~ send, and the recipient then has NI1U. • .UNIk = t~xi, i ]Ii=l slices to use within a reconstruction process. T h e procedure finishes when the recipient does not want more data, or all are sent. Matrix polynomials used in the sequel will be polynomials of the form A n x n + A ~ _ l x '~-1 + • • • + A 1x + A0, where x is a real scalar and all the coefficients Ai are members of C ~x ~, the space of complex r × s matrices. 2. M A T R I X INTERPOLATION,

POLYNOMIALS: DIVIDED DIFFERENCES

In this section, mathematical extension of Newton interpolation (fi'om scalar d a t a (xi, yi)) to a matrix generalization (for data (xi,Mi) with Mi E C rxs) is introduced. Let f : [a, b] --* C rxs be a matrix function and {x0 < xl < --. < Xn} be points belonging to an interval [a, b]. We will interpolate the function f and obtain a matrix polynomial of degree n, Pn(x), using the family of matrix polynomials { I , (X -- X O ) I , (X -- ZO)(X -- X i ) I . . . . , (X -- XO) " " " (X -- X n _ l ) I }

,

(1)

where the degree of the k th polynomial is exactly k - 1 and the leading coefficient is the matrix identity in C s x s, denoted by I. LEMMA 1. Let { K i ( x ) } ~ ° = o be a set ors x s matrix polynomials for which • K i ( x ) has degree exactly i, for i = O, 1, 2 , . . . , and • for each i, the leading coefficient of K i ( x ) is nonsingulax. Then, any r x s polynomial Q(x) with deg(Q(x)) _< n can be written uniquely in the form n

Q ( x ) -~ ~-~ C i K i ( x ) ,

/=0

where Ci ~ C rxs.

(2)

E.

306

DEFEZ et al.

PROOF. The following notation will be used:

Kl(x) = ho, Kl(X)

=

PllZ

Kz(x)

=

P2222

L(z)

= Pnnzn + pn,n-lz”-l +

&+1(z)

where ,&j E @““, deg(Q(z))

=

+ PlO, + P2ln: + P20,

Pn+l.n+l~n+l

‘..

+

+ L%L+l,n~n

+

&), ‘. . +

Pn+l,O,

0 I j < i, 0 5 i 5 n + 1, and ,& is nonsingular

= 0, then Q(X) = QO E FxS

and

Q(x) = (a~&,‘) Taking Cc = (ac&o1),

for 0 5 i < n + 1. If

Ko(x).

(2) holds in this trivial case.

If deg(Q(z)) = 1, then Q(Z) = orz + crc, with oc,ar and Cc satisfy the block linear system (w

Qo)=(G

Co)

E @“‘.

(?

Again, (2) is verified if Cr

ii:)7

but this has a unique solution because of the nonsingularity of prr and pee. Suppose, in general, that (2) is true for any polynomial of degree 1%2 1. A polynomial Q(X) = o,+r?+’

+ Q,?

+ . . + alx + ao,

cxi E C”“,

OLiln+l,

of degree n + 1, can be rewritten as follows: Q(X) = Q,+~x?+’

+ cry,9 +

‘. + ck!lZ + ck!o

= (-wn+K~,,,+,~n+1(4 + (-Q,+1P~~l,n+lP71+1,nZn -a,,+lP,-:l,,+lPn+l,o)

+ &zZn +. . . +

Kb,l(~)

= ~n+lP?&z+1

- (Y,+1P~~l,n+lP7L+1,n-121L-1 QlZ

- ”

+ a0

+ Rn(~),

where R,(z) is a matrix polynomial of degree at most n, and applying the induction hypothesis, we obtain that R,(z)

= 2 CzK(~), z=o

Ci E @“‘.

Thus, taking C,+r = ~~,+rPEir,~+r, (2) holds. Uniqueness of the representation (2) follows easily from repeated (TS times) use of a familiar +. . + Alz + A0 is identically the T x s zero, then scalar argument. If P(X) = A,? + An-lPdl each of its elements will be a polynomial in which all coefficients must be zero. That is. each T x s coefficient matrix Ai is zero.’

I

Let Pn(x) be the matrix interpolating polynomial [5] for the Crxs-valued function f(x) at the abscissae (50 < ~1 < ... < z%} c [a,b], and apply Lemma 1 using the family (1). Then, there exist matrices CO, Cl,. . . , C,, such that P,(X)

= c, + Cr(Z - IcrJ) + c&r

has the properties that Pn(Zi) = f(z,),

- Zc)(Z - X1) + . . . + C,(z for i = 0,1,2,.

- 20). ‘. (XT- x,-1)

. ,TL.

with deg(Ki(z)) = i, in which the leading ‘Lemma 1 does not generalize to the situation of a family {ITi(z coefficients, the flii above, are merely nonzero, rather than nonsingular. For example, with Ko(z) = [ “, i] and K~(x)

=

[i i]z, the nontrivial linear combination

[ : ~~]Kl(a) + [t I~]Ko(z) is identically zero.

307

Matrix Newton Interpolation

Hence, f (x0) = Co, f

(a)

=

f (2,) =

co

+ Cl

(51 -

20)

9

co+ Cl(GL - x0) +

With

(4 (x,-x0)... the linear system (3) can be written in block form as

(5)

where the coefficient matrix G is clearly nonsingular. Thus, the coefficients Cc,. . , C, are and vice versa. The next goal is to develop an algorithm uniquely determined by f(zc), . . . , f(q) for explicit determination of the coefficient matrices CO, Ci, . . . , C,. DEFINITION 2.

differences

Forf

recursively

: [u,b] --+ @“”

< 2,)

and {ZO < 51 < ...

C [a, b], define the matrix

fIXi = f(G), fixi,..

. ,

divided

by

xi+k~

=

fbi+l7.

i=O,l,..., . .7 2i+kl

-

xi+k

f[G, -

. .7 Zifk-I]

k ~ 1.

,

12,

(f-3)

xi

It is convenient, for future development, to write Ci = Ci[Xe, . . . , z,], i = 0, . . . , n, in order to emphasize their dependence on the values ~0, ~1,. . . , xn, so the interpolating polynomial, Pn(z), off at the points {zc,zi,...,z,} is

P,(x)

=~O[~O,...,~,]+~l[~Or...,~~l(~--o) +Cz[zo,...

) z,] (a: - 20)(2 - 21) + . . *

+ c, [x0,. . . 3Gzl (z - To). . . (x -

G-l),

with the Ci determined, one way, as forming the solution of the block linear system (5). We now seek a more efficient computational scheme for the Ci. First, consider a term Ck [Xi-i, . . . , &+k].

If D is any permutation of the indexes i - 1, i, . . . ,

i + k, G, [z&l,. holds. In fact, use the points {xi-i, mial P, Off at {Xi-i,. . . ,Zi+k}:

. . 7%+k]

=

. . . , zi+k}

ck

[zo(i-l),

. . . 7 xn(i+k)]

(7)

and construct the matrix interpolating polyno-

p(X)=C0[3Ji-1,...,Xi+k]fCl[Xi-l,...,Xi+k](X-X~,-1) +... fck[xi-_l,...

7%+k]

(x -

xi-l)

’ ’ ’ (x - %+k--1) .

E. DEFEZ et al.

308

constructin, c again, the matrix interpolating

Now, for the ordering {z,~~-.r~, . . . , ~,(~+k)}, nomial of f provides, from uniqueness,

P(x) = co[X,(,-l), . . . ,&(i+k)] + Cl + ”

[&T(,4)?~

+ ck

‘. &(i+k)]

[?+I),

. .

(x

1 “c(i+k)]

(x

-

poly-

Gr(z-1)) %(z-1))

”

(z

-

%(,+k--1))

.

Identifying corresponding terms with the preceding formulation yields (7). Next, let us verify that . ,xk]

ck-l[xl,. ck

for k = 1,2,3,.

(50,.

. . > zk]

=

-

. . , n. Construct the matrix interpolating

P(X) = co[Xo,.. .,Xk] +Cl[xO,..

.,xk](x

-x0)

On the other hand, with the points ordered as p(x)

ck-1

xk

= co[Xl,,...,Xo]

..rxk-l] (8)

polynomial P for f at {~a,.

ck[xO, .

+

+“’

[x0,.

x0

, zk](z

{zk, . . ,x0}, this same

+c1[5k,...,2O](LX--k)

-

20)

”

(x

. . . I zk]

=

(-20

+“‘+ck[zk,...,~O](~

-xl

-

.”

-zk-1)

. . ,x01(-x1-22-...-

ck[xk,.

=

ck

[x0,

. . . ,x1]

=

ck-l[zl,...,xk],

formula,

. ,zk__l]

+ Ck-l[Zk,.

x0]

ck[xk,...,

xk-1).

-xk)‘-‘(z--1.0).

+ Ck-_1[ZO,. xk)

-

P is

Identifying coefficients of xkP1, and applying a matrix version of the Cardano-Vieta Ck[xOr

, cck}:

.

(9)

,x1].

. . . >xk]

and Ck-l[xkr

hence, (9) can be written as ck[xO,.-,,

zk](-20

--xl

-“‘-~k-1)+Ck-~[50i...,Zk-l] =

ck[xO,...

fxk)

=

,xk]

(-Xl

-

22

-

.”

-

xk)

+ Ck-_1[Zl,.

,xk].

After simplification,

Ck[XO,*

‘.

,

xk]

(-x0

ck-l[xl,.

,Zk]

-

ck-1[ZO,.

.

,xk-I],

to yield (8). Finally, it will be shown that zk]

Ck[xO,...,

=

_f[xO,

.

(10)

, zk]r

for k = 0, 1,2,. __ , n. From the first equation of (3), Co[xo, . . .

)

z,] = f (X0)= flxol;

thus, from the second, G[~Or...r~nl

=

f

(21) Xl

-x0

f (20)

=

fbll

x1-ZO

fM

’

Matrix Newton Interpolation

and using

309

(6), f[xo

x1]

fhl - f[xol,

=

21 -x0 i.e., Cl[~O>...,

For the general

= f[zo,211.

step, if G-1

then,

%I

. .1 Zk-11,

[x0,. . . , Q-11 = f[zo,.

from (8), Ck[20,.

Ck-l[Q,. . . ,Ick]

. . ,~k]--k-l[~o,...,~k-ll 7

= xk

and taking

into account

the induction

so (10) indeed

[x0>

x0

hypothesis,

ck-l[zl,... ck-1

-

rxk] . .

, xk-l]

=

f[Zl,...rxk],

=

f[zO,

..

, xk-11,

holds.

3. The matrix interpolating polynomial P,(x) of the matrix function f : [a, b] -+ Vx”, at the points (x0 < x1 < . . ’ < 2,) c [a, b], can be written as

THEOREM

P,(x)

= f[xo] + f[zo,LQ](x

- x0) + ... + f[xo,.

. . ,GLl(x - 20) .‘. (x -

X?Pl)?

with coefficient matrices f[xo, . . . , Xk] given by (6). We refer to this as the matrix Newton divided difference form of the interpolating polynomial. Just as in the scalar case, the form for P,(x) appearing in Theorem 3 can be exploited construct an efficient algorithm, in essence, a matrix Horner’s rule, for evaluation of P,(x).

to

3. MATRIX POLYNOMIALS FOR PROGRESSIVE TRANSMISSION Section polation.

2 develops The scalar

Theorem

3.

a matrix situation

In this section,

polynomial generalization of Newton’s divided difference interindeed extends, to provide the interpolant in the form given by the matrix

polynomial

theory

is used to devise

an algorithm

for

progressive transmission of 3D images, and its visual effectiveness is illustrated in Figure 1. The matrix polynomial machinery thus gives an alternative approach for progressive transmission of images, to augment other schemes considered recently in this research and development area, matrix polynomials [5]. Addiincluding decorrelation techniques [3], wavelets [4], or orthogonal tional

approaches

may be found in [6,7].

Figure 1. A matrix polynomial interpolant, P,(z), progressive rendering. Newton interpolation at slices 41,43,45,. ,59, with Pg(z) evaluation at z = 42,44,46,. ,58. Single processor, PC time for interpolation and rendering: 3571.44 seconds.

310

E. DEFEZ et al

As discussed {(zi,M~)}~=a notation {($,

in the Introduction,

a digital 3D object

of n + 1, r x s matrices

assumed.

Now if pnk_i(z)

is frequently

stored as a set M =

Mi.

In the sequel, ~0 < ~1 < 22 < .. < 2, is the is the interpolating polynomial corresponding to nk pairs

Mf)}yLl which h ave been extracted from Ml, then evaluation of this intermediary PnI _i (x)

at all the n + 1 zis will provide a 3D reconstruction object Ml. Or, evaluation of P,l_i(z) for construction

of (i.e., an approximation

to) the original

at one value xj which is not one of the abscissae x” used

of Pnk_l (XT)creates a reconstruction

for an untransmitted

2D slice resident in

the original data base ML Within the overall process, appropriate ordering of all slices involved at any stage must be maintained, of course. Visual effectiveness of this process is illustrated below. The original data base, used in detail in Section 6, consists of 93 CT slices of a human head, each slice (matrix),

Mi, being of resolution 256

x

256. 2, = i for i = 1,2,3,.

. ,93 supplies a position

ordering (top to bottom of the head) for the corresponding slices Mt. The Mathematics

code is

given in Section 5, with the appropriate module in Section 5.2 for this part of the computations. To illustrate

direct use of the above matrix polynomial machinery, let F’s denote the inter-

polating polynomial corresponding to the transmitted (5% Mss). P 9( x ) was then evaluated at x = 42,44,46,.

data (41, Mdl), (43, Mds), (45, Mds), . , . ,58 to thus yield a combined stack of

18 slices for the skin rendering, given in Figure 1, of a head section which includes the nose. Even with this simple situation, processing time for the nine evaluations of P,(z)

using our PC system

described in Section 6 was approximately one hour. Even in scalar polynomial theory, repeated evaluations of higher-degreed polynomials can be CPU costly, or otherwise unsatisfactory, and one alternative to combat this is to use, instead, a piecewise-polynomial approximation consisting of lower-degreed segments. This is precisely the analog developed in the next section, to provide practicable PC-based computations in Mathematics for progressive transmission of 3D images using matrix polynomials. Two experiments are developed in some detail, and a piecewise linear structure is seen to be effective visually, as well as acceptable computationally.

4. PIECEWISE (LINEAR) MATRIX POLYNOMIALS FOR PROGRESSIVE TRANSMISSION A direct adoption of the preceding matrix polynomial machinery for progressive transmission, using Mathematics on a serial PC system, is not practical, as demonstrated in Section 3. Figure 1 therein illustrates that the machinery can, indeed, produce visually-effective progressive renderings, but at the expense of unacceptable elapsed times, at least using a single-processor PC system for the implementation. 2 A primary aim in this section is to extend the matrix polynomial development to piecewise polynomial structures and calculations. In [8], a matrix cubic spline approach produced substantial speed-up along with acceptable renderings. In fact, however, the simpler situation of linear segments has been found to still provide visuallyand , ,’ numerically-acceptable progressive renderings, and execution is faster. This is the development below: linear-piecewise matrix polynomial interpolants, with evaluations that exploit the matrix Newton divided differences, form the basis for practicable, single-processor, PC-based progressive transmission of 3D images. Full Mathematics code is given in Section 5. In [5], piecewise matrix polynomial interpolation, along with progressive transmission of subsets whose ordering was determined in a so-called balanced way, was used and it was concluded that desirable properties could be maintained in the rendered images even if segments consisting of just two subsets were joined together. Following notation in the Introduction, from a full image is selected for transmission (object) Ml = {(xi, Mi)}~zo, a subimage &lIk C M N (Ml U . . U i&-l) as the kth step in the progressive process, and Mi U. . .Ulw[k = { (zf , Mf)}yLl are the slices received by the recipient. 20ne extension

under consideration

is parallelization

of the computations.

Matrix

Thus, consider separation subsets as follows:

Newton

311

Interpolation

into components

of {(CC:, M,k)}y;r

which consist

of two-contiguous

jk, with j, E N such that jk + 1 = nk.

wherej=1,2,..., Now, apply

the Newton

but in each piece. Hence,

matrix interpolation, as was done in the above section (Theorem for Tt, j = 1,2,. . . , jk, Construct a hear matrix polynomial

P;(X) = f [xj”] + f [xj”,~j”+~](x

- x;)

j = 1,2 ,..., in the Newton

divided

difference

Mb1 -y "j+1 - xj

(x - XS),

(12)

j,,

form. Therefore,

Pk(x) = P!(x), is the piecewise

‘II; +

ZI

3),

matrix polynomial ,nk, and computing

if 5 E [~j”, zi+r 1,

step k. Thus, Pk(zk) = for the jth segment in transmission we obtain a full set of N slices, nk P”(Xi), i = 1,. . ,N, and n - n/, of them are generated by the matrix {M)}y$

MF, i =,l,... of them will be the original data polynomial P”(X). This will yield one basis for an approximate

rendering

image. As in the previous situation of use of a single polynomial, the recipient the ordering for all slices (received or computed) available, throughout

of the original

3D

software must maintain the entire progressive

process. 4.1.

Numerical

Often,

Errors

only visual

but we also include the progressive

assessment

is considered

some numerical

transmission,

for evaluating

measures

the goodness

using the Frobenius

define the overall

norm

of the reconstruction, 11. 1) [la].

At step k of

absolute error E, by

E = f: IlMj - P(xj)ll ,

(13)

j=o where P is the full interpolant. in Step k, denoted

Also, define the absolute error per pixel of progressive transmission

by EP, as EP=

r! rxsx(n+l)’

(14)

where r x s is the size of each of the n + 1 slices. We have found empirically, however, that E and EP are not useful measures for comparing results from differing methods, particularly if the data matrices, which represent 3D volume slices, have been constructed using different numerical scales. For example, VTK [9] uses an RAW format consisting of integers in the range [0,4095], whereas an RGB format has range [0,255] for each color component. Or, differing pixel resolutions might be used. So for future comparisons, we employ instead, the relative error of progressive transmission in step k, denoted by R, and the relative error per pixel of progressive transmission in step k, denoted by RP, as defined by

R=

nE JgoII4 II

and

(15)

312

E. DEFEZ et al R

RP= Of course,

if (zz, I&) is a slice sent, then

The full Mathematics numerical

/lMi - P(zi)ll = 0.

code for piecewise

error computations

is provided

(linear)

Auxiliary

polynomial

in Sections

5. MATHEMATICA 5.1.

(16)

r x s x (72+ 1)’ progressive

renderings

and relative

5.1-5.5.

IMPLEMENTATION

Functions

Fnorm computes the Frobenius norm of a matrix Fnorm [mat -1 : =Sqrt [Apply [Plus

, Flatten

[mat -21 I 1

EscribeRaw function writes a list datosto filenomFich in binary format EscribeRaw[nomFich_,datos_l:=BlockCstream, stream = OpenWriteBinary[nomFich]; WriteBinary[stream,datos,ByteConversion ->(ToBytes[#, ByteOrder->LeastSignificantByteFirst]&)]; CloseCstreamlI BU divides the list ep in sublists of 2 points where the last value of a sublist is the first in the next one BU[ep-]:=Block[qq=(Length[epl-I), Return [Table[ep [:I311, ep CCi+11111, i $1 ,qqlll busc&os

returns the position of the value J: in the subintervals determined

by the

list mat buscaPos[mat_,x_l:=Block[i=l, While [x>Last[mat[[illl&&icLengthCmatl,i++l; ReturnCill 5.2. Main Procedure Function for computing

the linear matrix Newton polynomial.

Oval is a list with 2

points and Lmat a list with two slices PolNewton[Lval_,Lmat_,x_l:=LmatCCIll+ ((LmatC ~21l-LmatC~1ll)/(Lval~~2ll-LvalCC11I>>* (x-LvalCClll) NW computes the Newton piecewise interpolating polynomial NW[x_l:=Block[i, i=buscaPos[values,x1; Return[PolNewton~valuesC~ill,imgsCCill,x11 5.3. Manipulating

I

and Preparing the Data

Loading package Utilities’BinaryFiles’.

EscribeRaw

inside functions need it

<
MatrixNewton Interpolation

313

5.4. Preparing the Experiment n=Length[images]; (*n=93*) (*slices transmitted in every step, for the set of *> (*93 slices with 256x256 pixels*) va1ues10={1,11,21,31,41,51,61,71,81,91,93}; va~ues30={1,4,7,11,14,17,21,24,27,31,34,37,41,44,47,51,54,57,61,64, 67,71,74,77,81,84,87,91,93}; va1ues50={1,3,4,7,9,11,13,14,17,19,21,23,24,27,29,31,33,34,37,39,41, 43,44,47,49,51,53,54,57,59,61,63,64,67,69,71,73,74,77,79, 81,83,84,87,89,91,93}; va1ues70={1,3,4,5,7,8,9,11,13,14,15,17,18,19,21,23,24,25,27,28,29,31, 33,34,35,37,38,39,41,43,44,45,47,48,49,51,53,54,55,57,58, 59,61,63,64,65,67,68,69,71,73,74,75,77,78,79,81,83,84,85,87, 88,89,91,93}; values90={1,2,3,4,5,6,7,8,9,11,12,13,14,15,16,17,18,19,21,22,23,24,25, 26,27,28,29,31,32,33,34,35,36,37,38,39,41,42,43,44,45,46,47, 48,49,51,52,53,54,55,56,57,58,59,61,62,63,64,65,66,67,68,69, 71,72,73,74,75,76,77,78,79,81,82,83,84,85,86,87,88,89,91,92,93}; (*Selection of the experiment. It must be changed in each series *> (*by values30, values50, etc.*) valuesi:=valueslO; Print[valuesil; Print["Percentage Print["Points

= " ,N[100*Length[valuesll/n]];

selected " ,valuesl];

(*Extracting selected images*) imgsi=TableCimagesCCvaluesl~~illll,i,i,Length[valuesl]]; Print["Points

selected ",valuesi];

(*Grouping the points*) values=BU[valuesl]; Print["Grouped

points ",values];

(*Grouping the images for interpolation*) imgs=BU[imgsl]; 5.5. Performing the Experiment (*Computation loop*) For[i=l;nor=O;ttot=O,i<=n,i++, xxx=Timing[NW[i]]; nn=Fnorm[xxx[[2]]-images[[il]] ; nor=nor+nn; ttot=ttot+xxxCCll1; Print["IMACE

",i,"

=== Norm error = ",N[nn,lO],

)' === Time = JJ,xxx[[i]]]; X2=Round[N[Flatten[xxx[C2]]1]]; (*Saving the interpolated images. *> (*Change 'NEWTON-IO' in each series*) EscribeRaw["NEWTON_lO.>'<>ToString[i] ,X2]

I;

314

E. DEFEZ et al

6. PC-BASED

EXPERIMENTS

AND RESULTS

The computations for the results described below were carried out on a PC Pentium III 450 Mhz with 128 MB RAM under Microsoft Windows 98 using Muthematica 3.0 1111,running the code given in the preceding section, for the piecewise-linear approach developed in Section 4. The 3D image was taken from the VTK CD-ROM

[9]. Of the many examples available, we

selected a set of 93 slices of a human head. There is a file for each slice, named from HEADSQ.l to HEADSQ.93.

The slices are B/W pictures in RAW format, with a B/W depth of 16 bits per

pixel and an image size of 256 x 256 pixels. The RAW format is not a format at all, but all the image information is stored, in this case, as an array of values between 0 and 65535, and it is necessary to know the kind of data being handled. On the other hand, this format is very easy for reading to and writing from Mathematics.

The 93 slices could be read using the Mathematics

sentence: ReadListBinary[NameFile,Intl6,65536, ByteOrder->LeastSignificantByteFirstl, where 65536 is the total size of the slice. After every slice was read and put in a list, the set of lists formed a matrix to be saved in a text file for further use. RAW format stores information pixel by pixel, one row after another (in this case, with a size 256 x 256 pixels) using 16 bits per pixel and unsigned values. Usually in CT images, the pixel depth is 12 bits, so here, the four most significant bits were not used, and the range of the entries of a CT image was therefore [0,4095]. On the other hand, it is interesting to note that the entries of a CT image do not represent color, but rather density, with increasing values for increasing densities. Once the calculations were done, we used a VTK program to render and display the 3D image created by our piecewise matrix polynomial interpolation procedure. 6.1. First Experiment

and Results

For the specific data set employed in the experiment and in terms of notation introduced at the beginning of Section 4, n = 92 and all data matrices M, are members of Q1256x256.5, = 0,. . ! 92 identifies the Mi, consecutively. This progressive transmission experiment consisted of sending sweepings of the 93 CT slices, in the following (balanced) way. 3 The first step, sending 11 slices, transmitted {(i, ~#Zc, 11/93 = 11.8% of the original information, an 88.2% compression rate, according to the selection scheme Mr: {1,11,21,31,41,51,61,71,81,91,93}. The second step, sending 18 slices more for a total of 29 slices (i.e., 31.2% of the original information or a 68.8% compression rate), yielded Ml1 u I& : {1,4,7,11,14,17,21,24,27,31,34,37,41,44,47,51,54,57, 61,64,67,71,74,77,81,84,87,91,93}. With the third step, transmitting 18 more slices, for a total of 47 slices (i.e., 50.5%) gave &!I1u &lI12 u M3 : {1,3,4,7,9,11,13,14,17,19,21,23,24,27,29,31,33,34,37,39,41, 43,44,47,49,51,53,54,57,59,61,63,64,67,69,71,73,74,77,79,81,83,

84,87,89,91,93},

3Selection for transmission among the 93 slices A4%can be automated in various ways. One being, for example, from following, in essence, binary subdivisions of the subscripts, with each sweep, as indicated by 46; 23,69; 12,35,58,81; 6,18,29,41,41,52,64,75,87; etc.

315

Matrix Newton Interpolation Table 1. Progressive 3D renderings, both skin and bone levels, using slices (resolution 256 x 256). First experiment. Captions show percent of the original data slices that have been transmitted.

Figure (a). 11.8%.

Figure (b). 31.2%.

Figure (c). 50.5%.

which is a 49.5% compression rate. For the fourth step, sending 18 slices more provided a total of 65 transmitted slices: MI

u R412U I&

U A4I4 : {1,3,4,5,7,8,9,11,13,14,15,17,18,19,21,23,24,25,27,

28,29,31,33,34,35,37,38,39,41,43,44,45,47,48,49,51,53,54,55,57, 61,63,64,65,67,68,69,71,73,74,75,77,78,79,81,83,84,85,87,88,89,

58,59, 91,93},

and this contained 69.9% of the original information, so a 30.1% compression rate had been

316

E. DEFEZ et al Table 1. (cont.)

Figure (d). 69.9%

Figure (e). 90.3%

Figure (f). Rendering of the complete original data set, 93 slices

achi evec . Finally, the fifth step, sending 19 slices more, to total 84 (i.e., a 90.3% data tram mi ssion, or a 9.7 I, compression rate) produced Mlr u Mz u MI3u &IL4 u Ms : {1,2,3,4,5,6,7,&g,

11,12,13,14,15,16,17,18,19,

21,22,23,24,25,26,27,28,29,31,32,33,34,35,36,37,38,39,41,42,43,

44,45>

46,47,48,49,51,52,53,54,55,56,57,58,59,61,62,63,64,65,66,67,68,

69,71,

72,73,74,75,7”,

77,78,79,81,82,83,84,85,86,87,88,89,91,92,93}.

A fter sending each set of subimages (Ml1, F&, IV&,Md, I&j), we computed a piecewi se- linear mat rix ,olynomial (P(s), P2(z), P3(z), P4(z), P5(2), respectively) using the matrix NIewton

317

Matrix Newton Interpolation divided puted

difference

approach.

for each j = 1,2,3,4,5

approximate

Then

a full manifestation

i = 0, 1, . . .92} could

be com-

to RAW format, for subsequent renderings of these For this renderin,, 0 the graphic modeler VTK [9,10] was

3D images of the original.

used, along with a program

{Pj(zi),

and translated

in TCL language

(131 that

can be found among

the examples

avail-

able with VTK, to allow one to compare, visually, the goodness and quality of a rendering. The renderings are presented at skin level and bone level in Table 1, for comparison with the original 3D object

renderings. Table 2. CPU time in seconds for computing and evaluating (12), to obtain interpolatory slices for the renderings in Table 1. First experiment. Compression

rate

CPU time

11.8%

31.2%

50.5%

69.9%

90.3%

128.22

84.9

52.87

29.53

9.44

The entries in Table 2 about CPU timing were obtained using the Mathematics function Timing[]. It is in order to note that in Table 2, we consider only the time for computing the interpolation calculations for unsent slices, i.e., 82 slices in Step 1, 64 in Step 2, 46 in Step 3, 28 in Step 4, and 9 in Step 5. Moreover, the time of the 3D rendering is not included at all, because the reconstruction

time is similar

in all cases, more or less 45 seconds,

Table 3. Numeric errors in the progressively transmitted images, vis-a-vis the original, at the five different data-use rates illustrated in Table 1. First experiment.

As illustrated above, visual assessment as transmission of CT slices progresses can be effective. Numerical computations are also informative for calculating errors in the decompressed image, as well as for summarizing and comparing processing times and storage requirements. Table 3 is self-explanatory: Section 6.2.

it gives numeric

values

of the error measurements

defined

by (13)-(16)

in

4.1.

Second

Experiment

and

Results

The original digital data form a block of 256 x 256 x 93 pixels, corresponding to a cuboid of actual size 384 x 384 x 279 millimeters. Hence, the resolution of the pixels is not the same along the three axes, X, Y, and 2. In fact, for the X and Y axes, there is a resolution of 1.5 millimeters per pixel, while for 2, it is three millimeters per pixel: 93 pixels (one pixel per slice) cover 279 mm along the axis 2, and 256 pixels must cover 384mm along each of X and Y. To explore differences, a second experiment was prepared by generating, with a simple C program, a new set of 256 slices of the head of resolution 93 x 256 pixels, from the old decomposition consisting of 93 slices, each of size 256 x 256, that was used for the first experiment. The 93 first slices of the head, see Section 6.1, can be visualized from top to bottom. These new 256 slices can be visualized from left to right. We repeated the experiment with approximately the same percentage of information sent, thus, more slices are needed in this second case to achieve a comparable bit level of information content. So, in this second experiment, n = 255 and all data matrices iMi are members of @g3x256. xi = 0, . . . ,255 identifies the Mi, consecutively.

318

E. DEFEZ et al.

In the first step we sent 27 slices, (11.5% of the 256 slices available) according to the selection scheme Ml : {1,11,21,31,41,51,61,71,81,91,101,111,121,131,141,151,161,171~ 181,191,201,211,221,231,241,251,256}. In the second step, sending 51 slices more for a total of 78, i.e., 30.5% of the data used, yielded Mi U ML2 : (1,4,7,11,14,17,21,24,27,31,34,37,41,44,47,51,54,57,61,64,67,71, 74,77,81,84,87,91,94,97,101,104,107,111,114,117,121,124,127,131,

134,137.

141,144,147,151,154,157,161,164,167,171,174,177,181,184,187,191,

194,197,

201,204,207,211,214,217,221,224,227,231,234,237,241,244,247,251,254,256}. For the third step, transmitting 51 more slices, for a total of 129 (i.e., 50.4% of the data), gave M~u~[ZU~~:{1,3,4,7,9,11,13,14,17,19,21,23,24,27,29,31,33,34,37,39,41,43,44,47, 49,51,53,54,57,59,61,63,64,67,69,71,73,74,77,79,81,83,84,87,89,

91,93,94,97,99,101,

103,104,107,109,111,113,114,117,119,121,123,124,127,129,131,133,

134,137,139,141,

143,144,147,149,151,153,154,157,159,161,163,164,167,169,171,173,

174,177,179,181,

183,184,187,189,191,193,194,197,199,201,203,204,207,209,211,213,

214,217,219,221,

223,224,227,229,231,233,234,237,239,241,243,244,247,249,251,253,

254,256}.

In the fourth step, sending 51 slices more yielded a total of 180 transmitted slices, 70.3%#,with IVlIiU Ml2 U IMaU &!I4: {1,3,4,5,7,8,9,11,13,14,15,17,18,19,21,23,24,25,27,28,29,31, 33,34,35,37,38,39,41,43,44,45,47,48,49,51,53,54,55,57,58,59,61, 69, 71,73,74,75,77,78,79,81,83,84,85,87,88,89,91,93,94,95,97,98,99,

63,64,65> 67,68, 101,103,104,

105,107,108,109,111,113,114,115,117,118,119,121,123,124,125,127,

128,129,131.

133,134,135,137,138,139,141,143,144,145,147,148,149,151,153,154,

155,157,158,

159,161,163,164,165,167,168,169,171,173,174,175,177,178,179,181,

183,184,185,

187,188,189,191,193,194,195,197,198,199,201,203,204,205,207, 214,215,217,218,219,221,223,224,225,227,228,229,231,233,234,235,

208,209,211,213, 237,238,239,

241,243,244,245,247,248,249,251,253,254,255,256}. And finally, in the fifth step, sending 51 slices provided a total of 231 transmitted slices? 90.2% of the original data, according to Mi U I& U Ms U M4 U Ms : {1,2,3,4,5,6,7,8,9,11,12,13,14,15,16,17,18,19,21,22,23, 24,25,26,27,28,29,31,32,33,34,35,36,37,38,39,41,42,43,44,45,46,

47,48,49,51,52,

53,54,55,56,57,58,59,61,62,63,64,65,66,67,68,69,71,72,73,74,75,

76,77,78,79,81,

82,83,84,85,86,87,88,89,91,92,93,94,95,96,97,98,99,101,102,103,

104,105,106,107,

108,109,111,112,113,114,115,116,117,118,119,121,122,123,124,125,

126,127.128,129,

131,132,133,134,135,136,137,138,139,141,142,143,144,145,146,147,

148,149,151,

152,153,154,155,156,157,158,159,161,162,163,164,165,166,167,168,

169,171,172.

173,174,175,176,177,178,179,181,182,183,184,185,186,187,188,189,

191,192,193,

194,195,196,197,198,199,201,202,203,204,205,206,207,208,209,211, 215,216,217,218,219,221,222,223,224,225,226,227,228,229, 236,237,238,239,241,242,243,244,245,246,247,248,249,251,252,253,

212,213,214, 231,232,233,234,235, 254,255,256).

Matrix Newton Interpolation

319

Table 4. Progressive 3D renderings, both skin and bone levels, using slices (resolution 93 x 256). Second experiment. Captions show percent of the original data slices that have been transmitted.

Figure (a). 10.5%.

Figure (b). 30.5%.

Figure (c). 50.4%

Then, following the same procedure as in the first experiment, but with these new data, the 3D renderings obtained appear in Table 4, with associated PC-computation requirements i in Table 5. Table 5 presents CPU times required for computation of unsent slices, for rendering-just 1% in the first experiment (for comparison), but here, 229 evaluations are needed for Step 1, 178 for Step 2, 127 for Step 3, 76 for Step 4, and 25 for Step 5. Again, the numeric values of the errror measurements, as defined in Section 4.1, are included, as Table 6.

320

E. DEFEZ et nl. Table 4. (cont.)

Figure (d). 70.3%.

Figure (e). 90.2%.

Figure (f). Rendering of the complete original data set, 256 slices. Second experiment. Table 5. CPU time (in seconds) for computing and evaluating (12), to obtain interpolatory slices for the renderings in Table 4. Second experiment. Compression CPU time

rate

10.5%

30.5%

50.4%

70.3%

90.2%

136.11

86.31

53.88

29.43

9.44

As in the error table for the first experiment, Table 6 gives numeric values of the error measurements defined in Section 4.1, but for this second experiment.

Matrix Newton Interpolation

321

Table 6. Numeric errors in the progressively transmitted images, vis-a-vis the original, at five different data-use rates illustrated in Table 4. Second experiment.

7. DISCUSSION The results

AND CONCLUDING

from the two experiments

described

in Section

REMARKS

6 provide

interesting

information

for comparisons and observations. A 3D digital object could be decomposed into a sequence of 2D subsets in many ways. For the test data set used here, the 3D head was decomposed into 93, 256 x 256-resolution parallel slices (top to bottom) for the first experiment and into 256, 93 x 256resolution parallel slices (left to right) for the second. The way in which slices will be selected for progressive transmission is, of course, important in the subsequent renderings. For the results reported here, only one scheme was used-repeated sweepings, in the increasing direction scalar ordering parameter, with a binary partitioning of this index within each sweep.4

for the

The way the object is decomposed into 2D subsets is also important, and this can be illustrated by the results from the two experiments. Visual difference is best seen by comparison between Figure (b) of Table 1 and Figure (b) of Table 4. For either the skin or the bone rendering, the latter pair (second experiment) simply “looks better”, even though approximately the same percentage of the original data was used in both cases, and the processing times were comparable. More specifically, for the first experiment’s pair, 31.2% of the original data was used with 84.9 seconds for processing. computing. 1.2362e-8 Table

In the second,

30.5% of the data

was incorporated

with

86.3 seconds

for the

The average relative error per pixel, however, was 1.6541e-8 for the first and only for the second, to provide a numerical measure, as well, for assessing Figure (b) of

1 versus Figure

(b) of Table 4.

This paper develops Newton divided difference interpolation for matrix polynomials (Section 2), for the progressive transmission of 3D images discussed in Section 3. The mathematical structure is adapted in Section 4, to provide a piecewise linear basis for progressive rendering of 3D images. With the adaptation, practicable Mathematics software is possible for implementation on a single processor, PC system, as timing Tables 2 and 5 show. The analysis and implementation presented here were based on a single 3D image, and 2D subsets of it. But, in fact, this restriction is unnecessary. Matrices in Sections 2 and 3 could be three dimensional (not two), so the object under consideration for transmission could be four dimensional (not three). This is an important extension currently for areas of computational medicine, for instance, functional magnetic resonance imaging captures 4D data as a succession of 3D images over time. Such MRI data could, of course, be “sliced” into a collection of subsets in differing ways. For medical diagnoses, 3D renderings at different instances in time constitute the most important way. For progressive data transmission, however, and perhaps for corresponding human interaction in focusing on subfeatures seemingly embedded within the data, other decompositions of the 4D object might prove more effective. The massiveness of such 4D data sets, combined with emerging needs for viable storage, compression, transmission and visualization algorithms, and software raise intriguing, open questions. We are currently investigating extension of matrix polynomial techniques for some of these. 4Current work underway addresses other schemes, including intelligent selection by the software itself. How to devise, a priori, a decomposition of the object that would provide a “good” rendering under progressive transmission of constituent subsets is a separate and intriguing question.

322

E. DEFEZ et al.

REFERENCES 1. S.G. Mallat, A theory for multiresolution signal decomposition: The wavelet representation, IEEE tins. PAMI 11 (i’), 84-95 (January 1980). 2. E.J. Stollnitz, T.D. DeRose and D.H. Salesin, Wavelets for computer graphics: A primer, Part I, IEEE Computer Graphics and Applications 15 (3), 76-84 (May 1995). 3. Y.-S. Kim and W.-Y. Kim, Reversible decorrelation method for progressive transmission of 3D medical image, IEEE mns. Medical Imaging 17 (3), 383-394 (1998). 4. E. Kofidis, N. Kolokotronis, A. Vassilarakou, S. Theodoridis and D. Cavouras, Wavelet-based medical image compression, Future Generation Computer Systems 15, 223-243 (1999). 5. E. Defez, A. HervL, A.G. Law, J. Villanueva-Oller and R.J. Villanueva, Progressive transmission of images: PC-based computations using orthogonal matrix polynomials, M&l. Comput. ModelZing 32 (lo), 1125-1140 (2000). 6. T. Sigitani, Y. Iiguni and H. Maeda, Progressive cross-section display of 3D medical images, Phys. Med. Biol. 44 (6), 1565-1577 (June 1999). 7. W. Wrazidlo, H.J. Brambs, W. Lederer, S. Schneider, B. Geiger and C. Fischer, An alternative method of three-dimensional reconstruction from two-dimensional CT and MR data sets, Med. Biol. Eng. Comput. 38 (2), 140-149 (March 2000). 8. E. Defez, A. Herv&, A.G. Law, J. Villanueva-Oiler and R.J. Villanueva, Matrix cubic splines for progressive transmission of images, Journal of Mathematical Imaging and Vision (submitted). 9. W. Schroeder and K. Martin, The VTK User’s Guide, Kitware, (1999). 10. W. Schroeder, K. Martin and B. Lorensen, The Visualization Toolkit. An Object-Oriented Approach to Graphics, Prentice-Hall, (1997). 11. S. Wolfram, The Muthematica Book, Wolfram Media, (1996). Computations, Johns Hopkins University Press, Baltimore, MD, 12. G.H. Golub and C.F. vanLoan, Matti (1989). 13. http://www.scriptics.com. 14. J.A. Ball, I. Gohberg and L. Rodman, Interpolation of Rational Matrix Functions, BirkhBuser, (1990). 15. E. Defez, L. J6dar, A. Law and E. Ponsoda, Three-term recurrences and matrix orthogonal polynomials, Utilitas Mathematica 57, 129-146 (May 2000). 16. R.C. Gonz&lez and R.E. Woods, Digital Image Processing, Addison-Wesley, New York, (1993). 17. L. JMar and E. Defez, Orthogonal matrix polynomials with respect to a conjugate bilinear matrix moment functional: Basic theory, Approx. Theory & its Appl. 13 (l), 66-79 (1997). 18. D.C. Kay and J.R. Levin, Graphics File Formats, Windcrest/McGraw-Hill, (1995). 19. M. Marcus and H. Mint, A Survey of Matrix Theory and Matrix Inequalities, Dover, New York, (1964). 20. K. Tzou, Progressive image transmission: A review and comparison of techniques, Opt. Eng. 26, 581-589 (July 1987).