Lecture Notes in Computer Science

A

(1)

Mode-1 Matrix

I

1

I

2

I

3

I

2

I

3

I

3

I

3

I

1

A

(2)

Mode-2 Matrix

I

1

I

2

I

3

I

3

I

1

I

1

I

1

I

2

A

(3)

Mode-3 Matrix

Fig. 1. Example of Unfolding the third order tensor to the three mode-n matrices

824

R. Xu and Y.-W. Chen

fold the third order to their mode-n matrices. The mode-n product of a ten-

sor

A ∈ R

I

1

×I

2

×...×I

N

and a matrix U

∈ R

J

n

×I

n

, denoted as

A ×

n

U, is an

(I

1

× I

2

× ... × I

n

−1

× J

n

× I

n+1

× ... × I

N

) tensor. Entries of the new ten-

sor is defined by (

A ×

n

U)

i

1

i

2

...i

n

−1

j

n

i

n+1

...i

N

def

=

i

n

a

i

1

i

2

...i

n

−1

i

n

i

n+1

...i

N

u

j

n

i

n

.

These entries can also be calculated by matrix product, B

(n)

= U

· A

(n)

,

B

(n)

is mode-n matrix of the tensor

B = A ×

n

U. The mode-n product has

two properties. One can be expressed by (

A ×

n

U)

×

m

V = (

A ×

m

V)

×

n

U =

A ×

n

U

×

m

V; and the other is (

A ×

n

U)

×

n

V =

A ×

n

(V

· U). The

inner product of two tensors

A, B ∈ R

I

1

×I

2

×...×I

N

is defined by

A, B

def

=

i

1

i

2

...

i

N

a

i

1

i

2

...i

N

b

i

1

i

2

...i

N

. The Frobenius-norm of a tensor

A is defined

by

A

def

=

A, A . The Frobenius-norm of a tensor can also be calculated

from its mode-n matrix,

A = A

(n)

=

tr(A

(n)

· A

T

(n)

), tr() is the trace of

a matrix.

3

Appearance Models Built by Generalized 3D-PCA

An medical volume, such as MR volume, can be seen as a third order tensor. The

three modes have their own physical meanings (coronal, sagittal and transversal

axis respectively). Supposing there are a series of medical volumes collected from

the same organs but different patients, there must be some common statistical

information of their appearance (the voxel’s intensity) for these volumes. In this

paper, we propose the multilinear algebra based algorithm named by G3D-PCA

to build such appearance models for 3D medical volumes.

We generalize the G3D-PCA for the appearance models building as follows.

Given a series of the third order tensors with zero-mean

4

A

i

∈ R

I

1

×I

2

×I

3

, i =

1, 2, ..., M ,

M

i=1

A

i

= 0, we hope to find three matrices with orthogonal columns

U

∈ R

I

1

×J

1

, V

∈ R

I

2

×J

2

, W

∈ R

I

3

×J

3

, (J

1

< I

1

, J

2

< I

2

, J

3

< I

3

) which is

able to reconstruct

A

i

from the corresponding dimension-reduced tensors

B

i

∈

R

J

1

×J

2

×J

3

with least error. Each reconstructed tensor ˆ

A

i

can be expressed by

Eq. 1. The matrices U, V, W contain the statistical intensity information of the

volume set

A

i

and they are respectively the orthogonal bases in the three mode

subspaces.

ˆ

A

i

=

B

i

×

1

U

×

2

V

×

3

W

(1)

We need to minimize the cost function S, shown by Eq. 2, to find the three

matrices.

S =

M

i=1

A

i

− ˆ

A

i

2

=

M

i=1

A

i

− B

i

×

1

U

×

2

V

×

3

W

2

(2)

4

If the tensors do not have zero-mean, we can subtract the mean-value from each

tensor to obtain a new series of tensors

A

i

which have zero-mean, i.e.

A

i

=

A

i

−

1

M

È

M

i=1

A

i

.

Appearance Models for Medical Volumes with Few Samples

825

Algorithm 1. Iteration Algorithm to Compute the Matrices U

opt

, V

opt

, W

opt

IN: a series of third order tensors,

A

i

∈ R

I

1

× I

2

× I

3

, i = 1, 2, ..., M .

OUT: Matrices U

opt

∈ R

I

1

×J

1

, V

opt

∈ R

I

2

×J

2

, W

opt

∈ R

I

3

×J

3

, (J

1

< I

1

, J

2

<

I

2

, J

3

< I

3

) with orthogonal column vectors.

1. Initial values: k = 0 and U

0

, V

0

, W

0

whose columns are determined as

the first J

1

, J

2

, J

3

leading eigenvectors of the matrices

È

M

i=1

(A

i(1)

· A

T

i(1)

),

È

M

i=1

(A

i(2)

· A

T

i(2)

),

È

M

i=1

(A

i(3)

· A

T

i(3)

) respectively.

2. Iterate until convergence

– Maximize S =

È

M

i=1

C

i

×

1

U

T

2

,

C

i

=

A

i

×

2

V

T

k

×

3

W

T

k

Solution: U whose columns are determined as the first J

1

leading eigenvec-

tors of

È

M

i=1

(C

i(1)

· C

T

i(1)

)

Set U

k+1

= U.

– Maximize S =

È

M

i=1

C

i

×

2

V

T

2

,

C

i

=

A

i

×

1

U

T

k+1

×

3

W

T

k

Solution: V whose columns are determined as the first J

2

leading eigenvec-

tors of

È

M

i=1

(C

i(2)

· C

T

i(2)

)

Set V

k+1

= V.

– Maximize S =

È

M

i=1

C

i

×

3

W

T

2

,

C

i

=

A

i

×

1

U

T

k+1

×

2

V

T

k+1

Solution: W whose columns are determined as the first J

3

leading eigenvec-

tors of

È

M

i=1

(C

i(3)

· C

T

i(3)

)

Set W

k+1

= W.

k = k + 1

3. Set U

opt

= U

k

, V

opt

= V

k

, W

opt

= W

k

.

In Eq. 2, only the tensors

A

i

are known. However, supposing the three ma-

trices are known, the answer of

B

i

to minimize Eq. 2 is merely the result of the

traditional linear least-square problem. Theorem 1 can be obtained.

Theorem 1. Given fixed matrices U, V, W, the tensors

B

i

that minimize the

cost function, Eq. 2, are given by

B

i

=

A

i

×

1

U

T

×

2

V

T

×

3

W

T

.

With the help of Theorem 1, Theorem 2 can be obtained.

Theorem 2. If the tensors

B

i

are chosen as

B

i

=

A

i

×

1

U

T

×

2

V

T

×

3

W

T

,

minimization of the cost function Eq. 2 is equal to maximization of the following

cost function, S =

M

i=1

A

i

×

1

U

T

×

2

V

T

×

3

W

T

2

.

There is no close-form solution to simultaneously resolve the matrices U, V, W

for the cost function S ; however if two of them are fixed, we can find the explicit

solution of the other matrix which is able to maximize S .

Lemma 1. Given the fixed matrices, V, W, if the columns of the matrix U are

selected as the first J

1

leading eigenvectors of the matrix

M

i=1

(C

i(1)

· C

T

i(1)

),

C

i(1)

is the mode-1 matrix of the tensor

C

i

=

A

i

×

2

V

T

×

3

W

T

, the cost function

S can be maximized.

826

R. Xu and Y.-W. Chen

I

1

I

2

I

3

I

1

J

1

I

2

J

2

J

1

J

2

J

3

opt

U

opt

V

opt

W

1

2

3

I

3

J

3

Fig. 2. Illustration of reconstructing a volume by the principal components

B

i

and the

three orthogonal bases of mode subspaces U

opt

, V

opt

, W

opt

Lemma 2. Given the fixed matrices, U, W, if the columns of the matrix V are

selected as the first J

2

leading eigenvectors of the matrix

M

i=1

(C

i(2)

· C

T

i(2)

),

C

i(2)

is the mode-2 matrix of the tensor

C

i

=

A

i

×

1

U

T

×

3

W

T

, the cost function

S can be maximized.

Lemma 3. Given the fixed matrices, U, V, if the columns of the matrix W are

selected as the first J

3

leading eigenvectors of the matrix

M

i=1

(C

i(3)

· C

T

i(3)

),

C

i(3)

is the mode-3 matrix of the tensor

C

i

=

A

i

×

1

U

T

×

2

V

T

, the cost function

S can be maximized.

According to Lemma 1, Lemma 2 and Lemma 3, we can use an iteration algo-

rithm to get the optimal matrices, U

opt

, V

opt

, W

opt

, which are able to maximize

the cost function S . This algorithm is summarized as Algorithm 1. In Algo-

rithm 1, we terminate the iteration when the values of the cost function is not

significantly changed in two consecutive times. Usually, the convergence is fast.

According to our experiences, two or three times iterations are enough.

Using the calculated matrices U

opt

, V

opt

, W

opt

, each of the volume

A

i

can be

approximated by ˆ

A

i

with least errors, where ˆ

A

i

=

B

i

×

1

U

opt

×

2

V

opt

×

3

W

opt

and

B

i

=

A

i

×

1

U

T

opt

×

2

V

T

opt

×

3

W

T

opt

. The approximation can be illustrated

by Fig. 2. In 3D PCA, the three matrices U

opt

, V

opt

, W

opt

construct the bases

on the third order tensor space; and the components of the tensor

B

i

are the

principal components.

Table 1. Comparison of G3D-PCA and the two PCA-based methods, given the training

volumes

A

i

, where

A

i

∈ R

I

1

×I

2

×I

3

, i = 1, 2, ..., M

Methods

Dimension of

Maximal Number of

Covariance Matrices

Available Bases

PCA–First Way (I

1

· I

2

· I

3

)

× (I

1

· I

2

· I

3

)

I

1

· I

2

· I

3

PCA–Second Way

M

× M

M

− 1

Mode-1 : I

1

× I

1

G3D-PCA

Mode-2 : I

2

× I

2

I

1

· I

2

· I

3

Mode-3 : I

3

× I

3

Appearance Models for Medical Volumes with Few Samples

827

Compared to the PCA-based methods, the G3D-PCA has advantages in two

aspects. Given M training volumes with the dimension of I

1

× I

2

× I

3

, the

covariance matrix for the first way to implement PCA has the a huge dimension

of (I

1

·I

2

·I

3

)

×(I

1

·I

2

·I

3

). So, this method is suffered from huge computing cost.

Different from the PCA-based method, G3D-PCA is to find three orthogonal

bases in the three mode subspaces rather than find one orthogonal bases in

the very long 1D vector space. The three bases are respectively obtained by

calculating the leading eigenvectors from three covariance matrices of the mode

subspaces. Since the three matrices have the dimension of I

1

× I

1

, I

2

× I

2

, I

3

× I

3

respectively, the calculation of eigenvectors is efficient. So G3D-PCA overcomes

the problem of huge computing cost. In the other aspect, the G3D-PCA can

obtain enough bases to represent the untrained volumes compared to the PCA-

based method implemented by the second way. In theory, the maximal number of

the tensor bases constructed from the bases of the three mode subspaces is I

1

·I

2

·

I

3

. So the number of the available bases is not limited on the number of training

samples and this make the G3D-PCA has good performance on generalization

even when few samples are used for training. Table 1 compares the differences

between the G3D-PCA and the two PCA-based methods.

4

Experiments

We use eighteen MR T1-weighted volumes of Vanderbilt [11] database to build

the appearance models of the human brain. These eighteen volumes are collected

from different patients, and their dimensions are 256

× 256 × 26. We choose one

volume as the template and align the other seventeen volumes onto the template

by similarity-transformation based registration. A 3D similarity-transformation

has seven parameters, three for translations, three for rotation angles and one

for scaling factor. Fig. 3 gives examples for three volumes of the registered data.

Fig. 3. Example of registered volumes

For these volumes, the first way to implement PCA is suffered from the huge

burden of computing cost. For the 256

× 256 × 26 volumes, the unfolded vectors

have the huge dimensions of 1703936, so the covariance matrix in the unfoled

vector space has an extremely huge dimension of 1703936

× 1703936. Suppos-

ing float type is used to store the covariance matrix in a computer, we need to

allocate the memory of 10816GB. This is impossible for the current desktop PCs.

828

R. Xu and Y.-W. Chen

The second way to implement PCA is suffered from the problem of bad perfor-

mance for generalization. This can be demonstrated by the leave-one-out exper-

iment, where the seventeen volumes are used to train the appearance model and

the left one is reconstructed by the trained models for testing. Fig. 4 gives the

result of leave-one-out experiment when the PCA is implemented in the second

way. It should be noticed that all the 16 available bases are used in the recon-

struction in this experiment. It can be seen that the quality of the reconstructed

slices is not satisfied. The reconstructed slices are blurred and the tumor part

on the 17-th can not be visualized from the reconstructed result. This problem

is because the available bases are too few to represent the untrained vector in

the unfolded vector space with large dimension.

7th-Slice

13th-Slice

17th-Slice

RMSE: 0.197

Original Slices

Reconstructed

Slices

Fig. 4. Result of leave-one-out testing for the PCA-based method implemented by the

second way

The G3D-PCA is not suffered from the two problems. The covariance matrices

for the three mode subspaces have the dimension of 256

×256, 256×256 and 26×

26 respectively, so the eigenvectors of them can be calculated efficiently. In the

other aspect, it is able to obtain enough bases to represent an untrained volumes

in the space of third order tensor by G3D-PCA. This can also be demonstrated by

the leave-one-out experiment. Fig. 5 gives the result of leave-one-out experiment

for the G3D-PCA. The training samples and testing volume are the same as

the experiment for PCA. The untrained volume is tried to be represented by

three tensor bases with different dimensions of 50

· 50 · 15 = 37500, 75 · 75 ·

20 = 112500, 100

· 100 · 20 = 200000. Since the number of the maximal available

bases is 256

· 256 · 26 = 1703936, the compressing rates for the three cases

are

50

·50·15

256

·256·26

≈ 2.2%,

75

·75·20

256

·256·26

≈ 6.6% and

100

·100·20

256

·256·26

≈ 11.7% respectively.

It can be seen that the quality of the reconstructed results become better and

Appearance Models for Medical Volumes with Few Samples

829

Original Slices

7th-Slice

13th-Slice

17th-Slice

Dimension of

the Bases:

15

50

50

RMSE: 7.69-E2

Dimension of

the Bases:

20

75

75

RMSE: 5.36-E2

Dimension of

the Bases:

20

100

100

RMSE: 4.23-E2

Fig. 5. Result of leave-one-out testing for the proposed G3D-PCA method

better with the increase of the dimensions of the bases. The root mean square

error (RMSE) is also calculated between the reconstructed volumes and the

original volumes. Compared to the PCA-based result, it can be seen that the

reconstructed results based on G3D-PCA is much better. The reconstructed

results are clearer. Additionally, the tumor region is also reconstructed well. This

experiment illustrates that the appearance models for medical volumes based on

G3D-PCA has good performance on generalization even when the models are

trained from few samples.

5

Conclusion

We propose a method named as G3D-PCA based on multilinear algebra to build

the appearance models for 3D medical volumes with large dimensions in this