Another type of measure is the percentage of topology violations among lattice

and voronoi neighborhoods.

Definition VII: Let Pneigh,l be the percentage of violations in the lattice

neighborhoods, and Pneigh,f, be the percentage of violations in the voronoi space

neighborhoods.

P = *100

neigh,l #(i:d (i,j)= l,i = ,..,K}

no

P = 100 (14)

neighxf #ti:d (i, j) = 1,i = 1,..,K

In cases where topology violation is inevitable, a SOM that contains minimum

Pneigh values is preferred. The last measure is the percentage of nodes participating in a

violation.

DefinitionVIII: Let Pnodes,l indicate the percentage of nodes that is a part of lattice

violation. Let Pnodes,f indicate the percentage of nodes that is a part of delaunay graph

violation.

sgn(H o,()) sgn(Ho ())

P = _*100 P = *-100 (15)

nodes,l K nodes,f K

Similar to Pneigh, a lower value of Pnodes is preferable, in cases where violations are

inevitable. By observing RD, Pneigh and Pnodes in several different SOM configurations, one

can make a quick decision about the suitability of a given SOM. Ideally, RD=I, Pneigh=0,

and Pnodes=0.

We performed tests on simple Id and 2d pattern spaces to demonstrate the use of

these three measures. The results are summarized in table 1. The first two pattern spaces

are the same, and consist of the contour of a circle. The succeeding two pattern spaces are

also same, but consist of not only the contour, but also the interior area of a circle.

Intuitively one would expect that the Id SOM mapping of a contour and 2d SOM

mapping of a full circle would be better than others. And, this is exactly what is

suggested by RD, Pneigh and Pnodes in table 1. The mapping of a full circle to a 2d SOM in

the last row indicate zero topology violation. On the other hand, the mapping of a circle

contour to a Id SOM in the first row demonstrate a slight violation. From the earlier

illustration in figure 12, it is known that this violation occurs only at the end points of the

SOM, because the circle is closed at these points, but the SOM lattice is not. This is

indicated indirectly by the low values of Pneigh and Pnodes. The values of Pneigh and Pnodes in

the middle rows are unacceptably high, indicating severe topology violations.

Given a SOM, the user can investigate the existence, nature, and acceptability of

topology violations using the measures presented here. In some cases, violations are

inevitable. However, choosing a SOM configuration with fewer violations is possible by

studying the proposed measures.

Table 1. Topology preservation quality of SOMs for ID and 2D data patterns

Data SOM ID RD Pneigh Pnodes Violation

size (%) (%) Type

Circle 1x16 1 1.87 6 12 Lattice

Contour

Circle 4x4 1 0.57 40 88 Lattice and

Contour Feature space

Circle 1x16 2 3.70 55 100 Lattice

Circle 4x4 2 1.0 0 0 None

There is an inherent shape-constraint imposed by the generic structure of the

SOMs. The generic shapes of the SOMs allow for open, non-intersecting curves and

surfaces. Standardization with respect to features that contain self-intersecting curves and

surfaces might cause topology violations, which in turn, jeopardize the alignment, and

create undesirable warps. Unfortunately, it is impossible to develop a cure for this

problem within the framework proposed in this study. It is recommended that this

methodology be applied exclusively to structures that are not self-intersecting, to avoid

topology violations completely.

CHAPTER 4

SELF-ORGANIZED ALIGNMENT AND WARPING WITH RADIAL BASIS

FUNCTIONS

The regularized alignment of individual features and warping the entire image is

addressed in this section. In our framework, SOMs are used exclusively for developing a

model for the feature during the alignment step. On the other hand, radial basis functions

are used to interpolate the new position assignments of individual image pixels during the

warping step. There is no coupling between the alignment and warping steps, so the SOM

and RBF interpolation procedures operate independently.

4.1. Self-organized Alignment

It is important to point out that, in the context of alignment of a pair of features,

regularization means imperfect matching of corresponding feature points for the sake of

smooth deformations in the surrounding area. There is a tradeoff between smoothness

and precision in the alignment. Regularization is accomplished at two levels in our

framework:

1) Vector quantization: Alignment is performed only at the anchor points obtained

through vector quantization. Precision is lost by not aligning all the feature points,

however, this loss in precision is compensated by a gain in smoothness. Examples that

demonstrate the benefits of vector quantization in terms of smoothness in interpolation

are provided in Bishop, 1995, chapter 5. As illustrated by Bishop (1995), the balance

between the loss in precision and gain in smoothness is achieved easily by adjusting the

number of anchor points.

2) Regularized matching of anchor points: At some anchor points, mismatches are

higher. If these anchor points are aligned as perfectly as others, the deformations in the

surrounding area will be larger than those in the other parts of the image. Therefore a

mechanism is needed to align the highly mismatched anchor points in a relaxed fashion.

This is the basic topic that is addressed in this section.

4.1.1. Background on Regularization Approaches

Regularized alignment of anchor points can be achieved through either

regularizing the SOMs or RBFs independently. Our literature search yielded only one

study (Goppert and Rosenstiel, 1996) which brings together regularization and SOMs. On

the contrary, regularization is a widely applied concept in interpolation with RBFs. In the

following, the existing regularization approaches in SOMs and RBFs will be summarized

briefly.

Regularized SOM lattice: At the end of SOM training, due to the final fine tuning

phase of SOM weight updates, in which the radius of the neighborhood reduces to one,

the positions of the codebook vectors become slightly disorganized, although their

ordering is kept intact. An extension to the SOM weight update rule can be made

(Goppert and Rosenstiel, 1996), so that the final positions of the nodes would reflect a

more regular lattice structure. Considering the SOM described in chapter 1, this extension

is formulated by replacing the weight update rule (2) with (16):

wk(t) = wk(t -1) + g(d, t)[xj wk( t -1)] + h(d,t) f(wk(t -1)) (16)

where g(d,t), and d are as described in chapter 1, f(.) is a multivariate function that returns

the ideal regularized position for Wk, and

SO, ifd = O

h(d, t) = _(2 / -2 with o set to either 2 or 3 to limit the range.

Xe elsee

Intuitively, this update rule works as follows. In the beginning of training, the

neighborhood function is much stronger, therefore the first term performs the

organization of the SOM. Towards the end, the effect of the first term diminishes, as

g(d,t) approaches zero. Therefore the second term becomes more important. At the last

stages of training, the winner is adapted towards the input vector, while its immediate

neighbors are adapted towards their regularized positions. In section 4.1.2 we adopt this

approach to retrain the template SOM in a regularized fashion for the feature in the

individual image.

RegularizedRBF interpolation: Given a set of sample points in the data space,

along with scalar target values, interpolation with basis functions is a process, which

associates scalar function values to all the points in the data space. This is an ill-defined

problem, because there are infinitely many solutions. In order to reach a pleasing unique

solution, smoothness constraints are introduced through regularization theory. The

smoothness criteria is enforced by the second term in the following error function, which

is minimized during interpolation.

K 2

E= L(d(v )-T(v)) + IP(T) 12, (17)

k=1 k k

where P is a differential operator, called stabilizer. The first term in the above equation

measures the accuracy of the function approximation, whereas the second term measures

the cost associated with the deviation from smoothness. The first term is small, when the

data, d, and the desired solution, T, are close. On the other hand, the norm of P(T) is

small, when T is a smooth warp. The regularization parameter, X, controls the

compromise between the degree of smoothness of the solution and its closeness to the

data (Poggio and Girosi, 1990a). Regularization theory restricts the types of functions

that can be used as basis functions to satisfy the smoothness criteria. Gaussian', and

inverse multiquadratic2 basis functions are supported by the regularization theory, while

multiquadratic3 functions are not (Haykin, 1999). The effect of the smoothness term in

the error function is reflected in the calculation of the coefficients of the radial basis

functions when ck are obtained by solving a set of linear equations. For the X translation,

the set of K linear equations (Poggio and Girosi, 1990b) are:

K -(01V -V 1 12/ya2)

T)x(vj)= Icck.(e k + S ) (18)

k=l k jk

= dx(vj) = vjx-wjx

where 67k is the Kronecker Delta (equal to 1 ifj = k; equal to 0 ifj k).

On the other hand, the smoothness term does not affect the calculation involved in

the interpolation of the whole image. The spatial displacements of each pixel, aj are still

computed in the x dimension by:

K -(lla -v 1/ k2)

Tx(aj)= Y cke

k=1

1 f(r) = exp(-r2/2o2), o>0

2 f(r) = 1 / (r2 + c2)1/2, c>

3 f(r) = (r2 + c2)1/2, C>

4.1.2. Proposed Method for Regularized Feature Alignment

In our framework, we have chosen to implement the regularization strategy in the

alignment step, through an extension to the SOMs. The reasons for this choice will be

clarified at the end of this section. There are two ways to develop a regularization

extension to SOMs. First, matching could be regularized with respect to the deformations

in the surrounding area of individual anchor points. Second, regularization could be based

exclusively on the positioning on the anchor points. In order to implement the first

approach, a region of interest, which consists of a subset of the whole image, must be

defined around each anchor point. At each iteration of the SOM, the part of the image

that falls inside the region of interest around the winning node must be warped and

evaluated for deformations. Accordingly, the SOM update term can be adjusted either

towards better fitting the individual feature, or for better smoothness. However, this

approach is very computational. It requires warping to be performed at every iteration,

which is quite time consuming. Therefore we developed the following regularization

scheme based on the displacements at the anchor points alone.

The basic error term at work in original SOM updates is the quantization error, Eq.

1 N 2

Eq= 111P -w |2 (19)

where 3n are the points on the individual feature, and w, is the closest node to 3n. The

quantization error measures how accurately the codebook vectors represent the feature

points. However, for regularized matching, accurate matching of the feature points is

unwanted at highly variable areas of the features. Instead, in these areas, less variability

between the atlas and individual features is preferred. Therefore an error term, Ed can be

constructed for measuring the discrepancy between the atlas and the individual feature.

1K 2

Ed= I v -w || (20)

Kk=l k k

Hence, the whole error term becomes:

E = Eq + Ed (21)

It is known (Kohonen, 1995) that the SOM algorithm does not employ a steepest

descent procedure to minimize E. However, it is shown (Kohonen, 1995) that the

resulting codebook vectors through the SOM procedure approximate the positioning that

would have resulted through a steepest descent procedure. Therefore, in order to

accommodate the new error term, we modified the original SOM weight update rule in

(2) as follows:

wk(t) = wk(t -1) + g(d, t)[xj wk( t -1)] + )(t)[wk(0) wk(t -1)], (22)

where wk(0) = vk(tend), as specified in the feature alignment section of chapter 1.The first

term in the weight update rule, (22), is the original SOM weight update term, which is

responsible for reducing Eq. On the other hand, the second term is for reducing Ed. The

parameter ) adjusts the balance between perfect matching of the individual feature, and

regularized matching. This approach is philosophically very similar to the approach taken

by Goppert and Rosenstiel (1996) which is described in the previous section.

It is important to note that there is a dilemma between accuracy versus regularity

(Goppert and Rosenstiel, 1996). Initially, when no action is taken, Ed is equal to zero,

because initially wk(0) = Vk(tend). As the SOM algorithm proceeds, Eq tends to be

minimized, while Ed increases. The final stage of the SOM training, in which the second

term dominates, aims to reduce Ed, and find a compromise between Eq and Ed. There are

several options to be considered in setting )(t). It could be a fixed small value, or a

slowly increasing function which reaches saturation towards the end of training. In our

simulations, we obtained successful results, particularly with the type of function

illustrated in figure 13. In figure 13, a(t) and )(t) obtain the same non-zero value towards

the end of training, and their values do not change from that point on.

t(t)

20 40

Figure 13. Sample training functions

The simulated 2d image warping example presented in figure 4 of chapter 1 is

repeated using the regularized SOM (RSOM) algorithm. The resulting individual SOMB,

and the warped image IB' are shown in figure 14a and 14b. As seen from figure 14b, the

deformations at the end points are more acceptable. The compromise for this smoothness

is made by reduced alignment accuracy. As seen in figure 14a, the individual RSOM

anchor points are positioned midway through the atlas and the individual SOM points.

Therefore, the loss in accuracy in RSOM can be predicted to be around 50% in

comparison to SOM alignment.

Using data from the above simulation, we investigated the change in Eq and Ed

during training. We also studied the distribution of error with respect to individual nodes,

after training. In figure 15a-15c, the change in Eq, Ed, and their sum are shown

respectively throughout training. As seen from these figures, the value of Eq in the

a) Location of nodes in b) Warped image using

SOM versus RSOM RSOM nodes

Figure 14. Alignment and warping results with regularized SOM

original SOM algorithm is lower than that of the RSOM algorithm. However, the reverse

is true for Ed, because the RSOM algorithm reduces Ed, while the original SOM

algorithm does not. Overall, the improvement in Ed is much higher than the improvement

in Eq, so the total error is less in the RSOM algorithm.

When the errors are plotted after training, at the individual node level, the figures

15d-15f are obtained. In these figures, individual nodes are indexed along the x axis via

their node id numbers. As seen from figure 15d, Eq is high for the RSOM at the end

points of the feature. This is due to the relaxed matching of the anchor points at the ends.

On the contrary, Ed is much higher for the original SOM algorithm, as seen in figure 15e.

Finally, as shown in figure 15f, the RSOM improves the total error twice as much

compared to SOM.

In addition to the above tests based on regularized SOMs, we have performed

tests with regularized RBFs, on simulated images. The performance of both

regularization approaches were similar, both computationally and in terms of smoothness.

Generally we did not find a performance improvement in one method over the other.

However, we noticed that regularized SOMs have two advantages:

1) The setting of the parameter, k is more intuitive, and easier, because it is in the

same domain with c(t), which is a part of the g(d,t) multiplier. The user needs to specify

the initial gain term, c(t), regardless of whether regularization is performed or not. This

facilitates the setting of the )(t) multiplier, which can be set in comparison to c(t), to

achieve a desired balance between accuracy and smoothness.

2) One basic assumption that we made is that the feature space is uniform. While

this is true for simulated data, it will not hold for real images, especially when the feature

is extracted automatically. SOMs are sensitive to the density of the data space, whereas

RBFs are not. Therefore for non-uniform features, the regularized SOM algorithm is

expected to perform better than the regularized RBFs, in obtaining a desired level of

smoothness.

One last issue to address is the complexity. The generation and retraining of the

standard SOM template both have the same complexity, O(n), for either regularized, or

non-regularized alignment, where n is the number of points in the feature space.

4.2. Proof of Preservation of Order Between the Template and Individual SOMs

The proposed alignment scheme is based on the assumption that the ordering of

the codebook vectors are the same in the template and individual SOM. If this

assumption does not hold, a valid correspondence will not be achieved between the atlas

and individual features. In this section it is shown that the order between the template and

individual SOMs is preserved, based on a condition imposed on the initial radius at the

RSOM

a) Quantization error, Eq, versus time

b) Deformation, Ed, versus time

d) Quantization error, Eq,

of individual nodes

e) Deformation, Ed,

of individual nodes

000 2 -C R

SOSOM

2-. SOM RO .- /1 RSOM

RSOM

c) Eq + Ed, versus time f)Eq + Ed of individual nodes

Figure 15. Errors Eq and Ed, during and after training

beginning of training. When the dimensionality is higher than id, due to the complex

interaction between the codebook vectors, the proof of SOM characteristics involve

intensive theoretical depth. Therefore we will outline the order preservation proofs only

for Id data.

a vi v2 VK b atlas

case 1c Wi W2 e WKd individual

iO--O wO'd

case 2 c wK WK-I e wi d individual

Figure 16. Id template SOM and probable orderings of individual SOM

Theorem 1: Let [a, b] be the Id atlas feature space, and {vi, v2, .., VK} be the

weights of the nodes of the template SOM as shown in figure 16. If a new SOM

containing weights, {wi, w2, ..., wK}, and initial radius, r, is trained using the individual

data set, [c, d], such that, wk(0) = vk, k=1,...,K and 1 < r < K, the ordering of all the nodes

Vk and wk are in topological correspondence.

Proof:

Part I, r < K: There are only two possibilities for the ordering of the individual

SOM, as illustrated in case 1, and case 2 of figure 16 (Kohonen, 1995). Case 1 is in

correct topological ordering with the atlas SOM, but case 2 exhibits reverse topological

ordering. We need to show that with the original SOM update rule, it is impossible to

obtain the ordering in case 2, when r < K. In order to prove that case 2 can not be

obtained, it is sufficient to show that after training, wi can not be positioned between

[e,d].

Consider feature points PE [e, d]. Since wK(O) = vK, and vK is the closest node to

points in [e, d], wK is the winner node for this set of P. Since r < K, only w2, w3, ...,wK are

updated with 3; wi never receives updates for PE [e,d]. Therefore according to Kohonen,

(1995), wi can only be positioned somewhere inside [c, e]. Indeed, if wi was the only

node to map into [c, e] it would have been placed at the centroid of [c, e]. But this is in

contradiction with case 2, which places wi outside [c, e]. Hence case 2 can not hold.

Part II, r > 1: Let r be 2. Initially, only a small subset of nodes will be winners,

because only few Vk are close to points in [c, d]. Assume only one node, wK, is a winner

initially. This will cause wK and wK-1 to be updated with P3e [c, d]. We know from

Kohonen (1995) that this will cause wK and wK-1 to be positioned somewhere inside [c,d].

But then this means that 3 3e [c, d], for which wK-1 is a winner, which in turn causes WK-2

to be updated with P. This way, not only wK, and wK-1, but also wK-2 will be positioned in

[c, d]. This process is repeated like a chain reaction, until all wk, k=l,...K are positioned

in [c, d].

Theorem 2: The regularized update rule in (22) does not spoil the ordering of the

codebook vectors, if )(t) is negligeably small initially.

Proof:

When )(t) is negligible, the original SOM update rule is in effect. It is known

from theorem 1, that in this case, Wk will be in the same order as vk. After this point, )(t)

increases. At this stage, the updates each node receives are multiples of vk-Wk, and the

first term in the weight update rule (22) is negligible. But this is equivalent to a reverse

mapping: from the individual, to the atlas (from wk e [c,d] to Vk e [a,b]). Since the

individual SOM's codebook vectors are in the same order with the atlas as shown in

theorem 1, this reverse mapping will initiate with the atlas ordering. Therefore, the

ordering of the nodes must remain the same.

4.3. Warping with Radial Basis Functions

Image warping is the final stage of the framework proposed in this study. Radial

basis functions centered at anchor points carry out the interpolation process for finding

the spatial displacements of each pixel in the atlas image. The equations involved in this

step are given in chapter 1. An important problem is determining the variances for the

Gaussian radial basis functions. If the variances are not chosen for optimal smoothness,

unacceptable warps, such as the ones shown in figure 18 are inevitable. The features in

figure 18 are the same as those presented in figure 4.

Haykin (1999) suggested an ad-hoc rule for determining the variances in an RBF

network, such that:

yk = dmax/ (2*K)1/2, k=l,..., K (23)

where dmax is the maximum distance between the anchor points, and K is the total number

of anchor points. However, this setting is sub-optimal and poor deformations may result.

For example, the warped image of figure 4 is generated by setting the variances as

suggested above.

4.3.1. Optimization of RBF Parameters

Ideally, the deformations in the warped image are expected to be local, but not

sharp. Small variances lead to sharp deformations, whereas large variances lead to global

deformations. If an energy term which can capture this feature of the deformations can be

found, then variances can be determined by minimizing the energy. We studied several

energy terms, including the energy of a spring, bending energy of a thin-plate, and the

energy of a rubber membrane. The energy of a rubber membrane was found to be the

most suitable, because its minimal values generated visually pleasing warps. Our

preference was reinforced after finding that some of the published nonlinear image

standardization procedures in the brain warping literature used the same type of energy

function (Ashbumer and Friston, 1999). The energy of a rubber membrane for a 2d image

is calculated using (24) below.

T (x,y) dT(x,y) dT (x ,y) dT (x ,y)

Em = x( 1 )2 +( )2"+( Y )2 +( Y y )2 (24)

m dox y ox ay

When we apply this energy definition to the interpolation term in (4), which is

used for assigning new spatial locations to pixels, we find:

K 2a (X -Xk + k 2)

E = I [- e k[( +y ( +

m k 21 k(

k 2b )2 2

2b -_( k 2)

+[- e y [(x +y)-(x +y)]] (25)

( k k

yk

where ak and 7xk are the coefficients, and variances for Tx, and bk and Gyk are the

coefficients, and variances for Ty.

Now we can study the suitability of this energy function. We would prefer the

value of the energy function to be high when the deformations are unacceptably sharp.

Also, as the deformations become more global, the energy term must increase. So, the

energy function should exhibit a graph with large values at low and high variances, with a

minimum somewhere in the middle. By setting all of the variances to the same value we

investigated the change in membrane energy on the simulated feature of figure 4. As

illustrated in figure 17, membrane energy has high values for very low and very high

variances. Optimum membrane energy is achieved only when the variances are allowed

30000

S25000

S20000

S15000

E 10000

a 6967=optimum

E 5000

0

0 10 20 30 40 50 60

variance

Figure 17. Change of membrane energy with respect to variance

to differ, and its value is half of that obtained through the suboptimal setting of the

variances. In figure 18, the effect of non-optimal setting of the variances is shown on the

simulated feature. When the variance is low, the deformation is unacceptably sharp, and

when the variance is high, global deformations are observed at the end points of the

feature.

a) Gk= 15 b) Gk =40

Figure 18. Two nonoptimal warps

We have used the steepest descent strategy to minimize the energy term in (25).

In figure 19a, 19c, 19b and 19d, results obtained from non-regularized alignment of the

linear features of figure 4 can be compared. Figures 19a and 19b are obtained using sub-

optimal variances, whereas figures 19c and 19d are obtained from optimal variances. The

warp in figure 19c is smoother than the one in figure 19a. The distribution of the

corresponding membrane energies in figures 19b and 19d illustrates that smoothness is

attained through a reduction in membrane energy. The total membrane energy Em is equal

to 11959.6 in figure 19b, and 6967.76 in figure 19d.

Similarly, in figure 19e, 19g, 19f and 19h, results obtained from regularized

alignment of the linear features of figure 4 are given. Figures 19e and 19f are obtained

using sub-optimal variances, whereas figures 19g and 19h are obtained from optimal

variances. The warp in figure 19g is smoother than the one in figure 19e. The total

membrane energy Em is equal to 2304.75 in figure 19f, and 658.45 in figure 19h.

Finding optimal variances for the warp is the most computational part of the

standardization procedure. Its order is O(P), where P is the number of pixels in the whole

image. In order to reduce the computations, the energy function might be evaluated not

on the whole image, but rather on a select subset of pixels (Woods et al, 1998a). The

pixels for which the evaluations are done can be spread randomly, uniformly, or more

closely around the feature. Assuming there are Q such pixels where Q<

computations can be reduced to the order O(Q), without significantly affecting the

outcome.

a) No regularization, Mj l warp

c) No regularization, optimal warp

e) Regularization, g m1 warp

d) E,, (total=6967.76)

f) Em (total=2304.75)

g) Regularization, optimal warp h) Em (total58.45)

Figure 19. Sample warps using SOM and RSOM for alignment

b) E, (total=11959.6)

4.3.2. Boundary Conditions

Up to this point, no restrictions are imposed for the new spatial assignments of the

pixels on the boundary of the image. In some applications, the outer boundaries of the

image, which consists of a rectangular frame for the 2d images, and a cubical box for the

3d images, must be immobilized. This can be achieved in three ways:

1) Penalty term: A new term can be introduced to the energy function to penalize

the situation in which the boundary moves. Therefore, the energy to be minimized can be

reconstructed as: E = y*Em + rl*Ep, where Em is the energy of the rubber membrane, Ep is

the penalty term, and y, fr are weights to adjust the effect of the penalty. In 2d, the penalty

term might be chosen as:

E, = Tx (xb yb2 y (xb,y )2 (26)

b

where b is a pixel on the boundary of the image.

2) Displacement multiplier: After finding the displacements of the pixels in the

whole image using Tx and Ty, these displacements can be multiplied with a 2d surface,

xy. Surface values are assigned to the image pixels such that, xy = 1 when the pixel

(x,y) is away from the boundaries. On the other hand, xy = 0, if(x,y) is on the boundary.

And when (x,y) is close to the boundaries, 0 < xy <1 holds. This way the boundaries are

pinned down with a smooth transition. We prefer to use this approach since it is the

simplest way to guarantee immobilized boundaries.

3) Boundary features: Boundary movement can also be avoided by defining them

as features. When the atlas image and the individual image have exactly the same

features, the discrepancies at the anchor points of these features is zero. Warping can be

accomplished by aligning multiple features rather than a single feature. Therefore, for 2d

images, 4 additional Id features can be defined, consisting of the sides of the bounding

frame. For 3d images there must be 12 additional 2d features, one for each side of the

bounding cubical box. The treatment of multiple features, as summarized in the next

section, is merely a warping issue, which does not affect the alignment step.

4.3.3. Generalization to Multiple Features

In this section, we show how to handle the standardization problem with respect

to multiple features. First, we need to revise the definition of features, based on the

definition of a single feature given in chapter 1.

Definition I:

Let IA and IB be the atlas and individual images respectively. Consider two sets of

structural features, AE IA, and Be IB in these images, consisting of U elements, such that:

A = {Ai, A2, ...., Au}, and B = {B1, B2, ...., Bu},

where A1 = {y1, U12, ..., (} and B1 = {3n, 312, ..., I1N},

A2 = {021, U22, ..., 2p} and B2 = {21, 322, ..., 32Q}, ... etc.

Here, tim, m = 1,.., M and Pin, n = 1,.., N are points in 2d or 3d space in which IA

and IB are defined, and M and N are not constrained to be equal. The same argument

holds for the points aum, and (Xun in the other features Au and Bu.

The objective is to align Au with Bu, u =1, ..., U by generating anchor points, and

to warp the image IB by interpolating the discrepancies between the whole set of anchor

points. The alignment step is not affected by having multiple features, because alignment

is a process which is exclusive to a pair of features, one from the atlas image, the other

from the individual image. However, we still need to revise the notation of the data

structures in the template generation and alignment steps.

Definition II:

Let the anchor points of atlas feature, Au, be Vu = {vul, Vu2, ..., VuK} and the

anchor points of the same feature, Bu, in the individual image be Wu = {wul, Wu2, ...,

WuK}. Define the discrepancies as Du = {vul ul, Vu2 Wu2, ...., VuK WuK}.

The multiplicity of the features only affect the warping step. Since warping is

accomplished by interpolating the discrepancies at the anchor points, there is no reason to

discriminate between the anchor points as belonging to specific features. The set of

displacements, Du, u = 1, ..., U can be pooled under one displacement set, D. Therefore,

the RBF coefficients must be calculated for this whole set, and variances of the RBFs at

all anchor points must be optimized at the same time. The spatial displacements for

individual pixels can be found as follows.

Let T be a 3D transformation: T = (Tx, Ty, Tz), where Tx(x, y, z), Ty(x, y, z), Tz(x,

y, z), are real valued multivariate functions, R3 --R. Compute d(aj), the displacement

vector, at voxel aj e IA, by:

d(aj) = (Tx(aj), Ty(aj), Tz(aj)), where

K -(|a -v H12/o 2) K -(lla -v 112 /a 2)

T (a)= clk.e k k + C22.e C2k 2k

x k=l k=l

K -(a -v l12/y 2)

+ IUk.e J (27)

k=1

Here Cuk, and Guk are the coefficients and the variances for the RBFs in feature Au

respectively. Ku is the number of anchor points. The coefficients can be obtained by using

the known displacement vectors, duk, at the anchor points, solving the set of K + K2 +.

+ Ku linear equations:

Tx(Vuk) = (Vuk-Wuk)x4 (28)

4 (.)x is an operation that returns the x component of the vector in the argument

CHAPTER 5

3D APPLICATIONS: STANDARDIZATION OF MRI OF THE BRAIN

One important application of neuroimaging is localization of brain activity. In

localization studies, neuroanatomical correlates of function are investigated by overlaying

the functional activity in multiple brains. However, the anatomical variability hinders

precise localization of the activity, as illustrated previously in figure 5 of chapter 2.

Reduction of normal variability non-linearly, without introducing sharp deformations is

the ultimate objective in standardization of brain images.

While the current automatic methods are quite efficient in normalizing the overall

size and shape of brain, mismatches of up to 1.5 cm in major anatomic structures such as

the sylvian fissure, central sulcus, parieto-occipital sulcus and superior temporal sulcus

still prevail after standardization (Grachev et al, 1999; Steinmetz et al.; 1989, Leonard et

al., 1998). Considering the complexity of the anatomy and the extent of variability in

individual structures, this result is not surprising. Although most structures follow a

general pattern of curvature, they also exhibit broad differences in continuity, branching,

and positioning. Further details of the types of variability in major cortical structures can

be found in Ono (1990) and Royackkers (1999), and numerical figures are provided in

Rademacher et al. (1993), Steinmetz et al. (1989), and Thompson et al. (1996a).

In the earlier chapters, we have presented a framework for performing smooth

non-linear alignment of the individual features. However, there are two important issues

to be resolved before applying these procedures to brain standardization.

1) Selection of the standard image:

Currently, there is no consensus in the neuroimaging field regarding the standard

atlas brain. The widely accepted Talairach atlas (Talairach and Tournoux, 1988) consists

of manual sketches of the brain of a 57 year old woman. Although the coordinate system

of the Talairach atlas is inarguably the current standard reference system, the validity of

the atlas itself as a standard is not yet established. This is because the atlas does not

consist of MR images, and a single brain can not exhibit the average characteristics of all

brains. If the chosen atlas brain contains average characteristics of the structures, at least

it is guaranteed that the atlas itself does not contain some extreme features. In addition,

the discrepancies between the atlas and any given brain will be less pronounced.

Therefore, a widely accepted approach is the creation of probabilistic atlases derived

from multiple brain images. In the past six years, following an initiative for the creation

of probabilistic brain atlases (Roland and Zilles, 1994; Mazziotta et al., 1995), several

groups independently generated probabilistic atlases consisting of merged topologies

from 10-300 brains (Evans at al., 1994; Roland et al., 1994; Thurfjell et al., 1995; Woods

et al., 1999; Thompson et al., 2000).

On the other hand, in the literature, it is a general practice to test new

standardization methods on a standard image which is randomly chosen from the image

pool of that particular study. Obviously, such an image cannot be accepted as a general

standard atlas. However, this choice facilitates quick and easy use of the newly developed

methodology, and helps demonstrate the benefits of the proposed method without getting

into the format and access difficulties of a third-party atlas image. In the standardization

procedures in this chapter, we subscribed to this approach, and dedicated one randomly

chosen brain MRI as the standard image. We are aware that such a choice of the atlas is

not ideal, but it is sufficient to demonstrate the results of our method. In the future, the

studies presented in this chapter can be replicated on a larger subject pool, and a

generally accepted atlas image.

2) Quantitative evaluation of the variability:

Evaluation of the variability in brain images is another subject on which

consensus is not yet established partly because the optimal metric depends on the

particular application. For an efficient evaluation and benchmarking of the

standardization algorithms, differences in the four following issues need to be resolved.

a) Definition of error: Several different types of errors are suggested in the

literature. Broadly these can be grouped under two categories: registration errors and

volumetric errors. Registration error is defined as the distance between an anatomic

landmark on the standard image and on the warped individual image. Ideally, after

standardization, the landmarks on the standard and the warped image should coincide,

resulting in zero registration error. Registration errors may differ according to the

distance term used. On the other hand, volumetric errors can be used to capture the

variability more generally, usually along an anatomic structure, rather than distinct

anatomic landmark points. As illustrated in figure 20, three types of mismatches are

possible while comparing the volumes of the structure in the standard and the warped

image. Type I: The feature on the standard image has points that lie outside the feature on

the warped image. Type II: The feature on the warped image has points that lie outside

the feature on the standard image. Type III: The features on the standard and warped

image are completely mismatched.

Type I Type II Type III

Figure 20. Volumetric mismatch types

For instance, Gee et al. (1993) defined volumetric discrepancy as:

V r)V

A B (29)

V uV

A B

where VA is the volume of the given structure in the standard image and VB is the volume

of the same structure in the warped individual image. Ideally, when the structures in the

standard, and the warped image completely overlap, this volume difference becomes one.

While it captures both type I and type II of mismatches, the volume error definition in

(29) does not indicate the severeness of type III mismatch.

Another volume definition by Fischl et al. (1999a) is as follows:

V -V

B,all Bavg (30)

V

B,avg

where VB,all is the union of the volume of the given structure across all subjects and VB,avg

is the average volume of the given structure among all subjects in the warped images.

Ideally, if the volumes overlap completely, this volume error will become zero. Only type

I and type II errors are captured in this formula for multiple subjects.

In order to report volumetric errors, we preferred to define a new volume error

which has ability to capture all three types of mismatches. The details of this volumetric

error are provided in section 5.3.

b) Selection of landmarks: In Grachev et al. (1999), an extensive set of anatomic

landmark points (128 per hemisphere) is defined. Several other sets of landmarks are

used throughout the literature to demonstrate the results of individual standardization

algorithms. However, there is no study that brings together all reported landmarks to

validate their use and summarize the underlying incentives in their selection. While

reporting the registration errors, we selected landmark points that are relevant to the

structures in our study using the list of landmarks provided in Grachev et al. (1999).

c) Measure of deformation: The distribution and scale of the deformation that the

individual images undergo during the warping process are very important measures.

Some of the standardization algorithms in the literature, including ours, take the energy

embodied in the deformation into account while choosing the parameters of the warping.

However, there is no numeric report in the literature about the scale of deformation, nor

its change relative to the parameters that are assigned independently. In the context of

brain images, it is especially important that the deformation is smooth and minimal. At

the end of section 5.3, we provide the scaling of the deformation in our framework with

respect to the only independent parameter, the number of SOM nodes.

d) Standard test set: In order to compare the performance of different

standardization techniques reliably, all algorithms must be run on the same data set

consisting of the standard image and multiple individual images. Currently such a test set

is not established for brain imaging applications. The research reported in the literature

are all performed on independently chosen test sets. In addition, the scan parameters also

vary from one test set to the other. Therefore, we have collected our own test data as

explained in the next section.

5.1. Data Collection and Preprocessing

Multi-spectral MRIs consisting of TI and T2 contrasts are collected from 10

subjects. Our motivation behind this type of image acquitision was to enable the

preprocessing of the same data set through both automatic and manual feature extraction.

In the literature, multi-spectral MRIs are shown to be efficient in automatic segmentation

of CSF, gray matter and white matter (Hall et al., 1992; Bezdek et al., 1993; Clarke et al.,

1993; Fletcher et al., 1993; Ozkan et al., 1993; Lundervold and Storvik, 1995; Alfano et

al., 1997; Vaidyanathan et al., 1997). Once these three types of tissues are segmented,

specific features can be extracted using shape recognition techniques.

The MRI acquisitions are performed on a 1.5T Siemens with the following scan

parameters. For the T1 weighted images, TR=15 ms, TE=7ms, FA=150, slab number=l,

slab thickness= 160mm, number of partitions=160, shift mean=0mm, imaging

plane=sagittal, field ofview=250mm, voxel size=0.98x0.98x1.00mm, scan time=10 min.

For the T2 weighted images, TR=3300ms, TE=105ms, FA=1800, slab number=8, slab

thickness=20mm, number of partitions=16, shift mean=0mm, imaging plane=sagittal,

field ofview=250mm, voxel size=0.95x0.98x1.25mm, scan time=19.5 min.

5.1.1. Transformation into Standard Coordinates

An affine transformation is performed in order to transform the images into

Talairach coordinates, which is considered to be the current standard coordinate system

for the brain. In order to perform the transformation, 10 anatomic landmark points1 are

specified manually using the AFNI toolbox (Cox, 1996). Using these landmarks, the

1 List of landmark points: Anterior commissure, posterior commissure, two points on the interhemispheric

fissure, two extreme points on the brain for each dimension (a total of six points)

affine transformation brings two structures, anterior commissure (AC) and posterior

commissure (PC) into register with the standard coordinates, aligns the line between AC

and PC with the horizontal axis, and coarsely scales the brain using the specified 6

extreme points on the outer brain. In figure 21, a brain image is shown before and after

the affine transformation into Talairach coordinates.

..': :. : .' . . .... .....: .. . .

t 1..*

Figure 21. A sample brain image before (left) and after (right)

transformation into standard coordinates

5.1.2. Feature Extraction

All the features presented in this dissertation, including some of the simulated

ones are extracted manually using LOFA, a software developed in our lab (Gokcay,

1999). LOFA, Localization of Functional Activity, is a toolbox written in Pv-Wave for

delineating anatomic structures in individual subjects, and extracting the corresponding

region of interests in overlaid functional images. Additional information on the

operational principles of LOFA is provided in the Appendix.

The benefits of self-organizing feature based alignment is maximized in the

standardization of brain images, if several structures are chosen to represent the entire

brain. Some of these structures are Sylvian fissure, superio-temporal sulcus, parieto-

occipital sulcus, cingulate sulcus, central sulcus. Unfortunately, there aren't many sulci in

the frontal areas of the brain that follow a general pattern of curvature having variability

in the ranges comparable to the sulci listed above. The variability in the frontal lobe is

higher than in other regions, and therefore the deformations in this area would be greater,

if frontal structures are included in the standardization.

In this section, our main objective is to compare the performance differences of

the proposed self-organizing framework with respect to the alignment of topologically

different structures. Therefore we have chosen two structures, the Sylvian fissure and the

cingulate sulcus to demonstrate our method. The Sylvian fissure starts on the lateral

aspect of the brain and runs deep. Most of the current standardization algorithms result in

high registration errors for landmarks on the Sylvian fissure. On the contrary, the

cingulate sulcus is a thin structure located on the medial aspect of the brain. According to

the figures in Steinmetz et al. (1989), the maximum variation in the posterior part of the

Sylvian fissure and the marginal ramus of the cingulate sulcus are both 2 cm. However,

when variability is compared in all areas of these two sulci, the cingulate sulcus is aligned

somewhat better during transformation into Talairach coordinates, since this

transformation specifically aligns two other midline structures AC and PC.

While tracing these two structures using LOFA, we established some guidelines

for consistency. The cingulate sulcus is traced sagittally, slice by slice, starting in the

midline, working laterally, until it no longer is visible. On the average, there were 13

slices that qualified. Whenever the sulcus is interrupted, the traces are discontinued and a

new trace is started at a new branch. Traces are specifically positioned on the sulcus,

rather than the gyral area. The posteriormost point is chosen as the marginal ramus, and

the anteriormost point is chosen as the point that meets a vertical line passing through the

genu of the corpus callosum. Sample tracings of the right cingulate sulcus (RCS) are

shown in figure 22a.

Similarly, the Sylvian fissure is also traced on the sagittal slices. The drawings

proceeded from medial to lateral. The first slice is chosen as the one where the Heschl's

gyms is clearly visible, and the insula appears in the background. The last slice is

determined as the lateralmost slice where the frontal brain disappears at the background

and only the temporal lobe becomes visible. On average, there were 24 slices that

qualified. Traces are specifically positioned on the lower bank of the Sylvian fissure,

other branches were left out. The posterior and anteriormost points are found at the two

end points of the horizontal ramus. Sample tracings of the right sylvian fissure (RSF) are

shown in figure 22b.

5.2. Standardization

As summarized in chapter 1, standardization is performed in two phases. In phase

I, the SOM templates are generated from the structures in the atlas image. In phase II,

these SOM templates are retrained using the structures in the given image for alignment,

and warping is performed using the difference vectors generated during alignment. In this

section, the steps involved in both phase I and II are discussed.

5.2.1. Template Generation

In order to choose the correct dimensionality for the SOM templates of the two

features, RCS and RSF, the topology evaluation computations presented in chapter 3 are

performed. Results of these evaluations are shown in table 2. As seen from this table, the

correct dimensionality for RCS is 1, and RSF is 2. Because otherwise, severe topology

viloations are encountered, as indicated by the Pnodes value of the rows two and three. The

slight topology violation observed in the 2d RSF is investigated. This violation is

consistent across the representation of Sylvian fissures from different subjects, hence

does not present a conflict in the one-to-one correspondence. The violation occurs

because of the flattening of the Heschl's gyms as the slices move from medial to lateral.

". ". :" ...

a) RCS, first eight medial to lateral slices ordered first left-right, then top-bottom

b) RSF, first eight medial to lateral slices ordered left-right then top-bottom

Figure 22. Tracings of two sulci (RCS, RSF)on the standard image

Table 2. Evaluation of topology preservation of the RCS and RSF SOM templates

Data SOM ID RD Pneigh Pnodes Violation

size (%) (%) Type

RCS 1x10 1 1.0 0 0 None

RCS 7x3 1 0.5 34 76 Feature space

RSF 1x20 2 2.1 34 75 Lattice

RSF 7x3 2 0.85 6 19 Feature space

a) Id SOM model (1x20) of RSF

c) Id SOM model (1x10) of RCS

b) 2d SOM model (3x7) of RSF

d) 2d SOM model (2x10) of RCS

Figure 23. SOM templates for two sulci: Right sylvian fissure (RSF)

and Right cingulate sulcus (RCS)

In figure 23, the SOM templates that are evaluated in table 2 are visualized. As

seen in figures 23b and 23c, the 2d RSF and Id RCS result in perfect lattices with equally

spaced nodes in the feature space. On the other hand, the Id RSF in figure 23a, and the

2d RCS in figure 23d exhibit nodes wrapping around in the feature space and vice versa.

5.2.2. Self-organized Alignment and Warping

Alignment of the features in a given brain image is achieved by retraining the

SOM templates of figure 23b and 23c with the points belonging to the individual

structures of the given brain image. During retraining, we used an initial radius of 2, and

the initial gain and regularization multipliers are set as 0.2 and 0.05. In figure 24, the

difference vectors are visualized for each SOM node after alignment. As seen in figures

24b and 24d, regularization improves the uniformity of the discrepancy vectors.

In figure 25, the resulting standardized images are shown for one slice in each

structure, RCS and RSF. The tracing of the atlas structure is overlaid on all images to

indicate the differences between the feature on the standard, given and warped images.

As seen in figure 25b, there are significant differences between the atlas and the given

images in terms of the positioning of the two features, RCS and RSF. After warping, the

differences are reduced as seen in figures 25c and 25d. The matching of the features after

SOM alignment is better in comparsion to RSOM alignment. However, SOM alignment

results in high deformation at the anterior tip of the Sylvian fissure, whereas the

deformation in RSOM alignment is smoother, hence, more acceptable. Numeric

comparison of membrane energy and matching errors for the standardization of all 10

brains is provided in section 5.3.

84

Finally, in figure 26, the features from all 10 brains are overlaid on the standard

image to visualize the discrepancy after transformation into standard coordinates and

warping simultaneously. As seen in figure 26a, the registration error in some parts extend

to 1.5-2.0 cm. This error roughly drops by half after applying the framework proposed in

our study. SOM alignment provides better matching as indicated by the difference in the

alignment of figures 26b and 26c. However, the advantage of RSOM alignment in

smoothness is not clear. This is covered in figure 27c of section 5.3.

a) Discrepancy in SOM alignment of RSF

lBo-

1^0 -

: ^

b) Discrepancy in RSOM alignment of RSF

c) Discrepancy in SOM alignment ofRCS d) Discrepancy in RSOM alignment of RCS

Figure 24: Discrepancy vectors for two sulci:Right sylvian fissure (RSF)

and Right cingulate sulcus (RCS)

i

B

i

c

t

~ i-

a) Standard brain (IA)

b) Individual brain (IB)

c) Warped individual brain (IB') with RSOM alignment

d) Warped individual brain (IB') with SOM alignment

Figure 25: Sample warps with SOM and RSOM alignment on the

RCS(left) and RSF(right)

a) After affine alignment

b) After RSOM alignment

c) After SOM alignment

Figure 26: RCS (left) and RSF (right) traces overlaid from 10 brains after

affine, RSOM, and SOM alignment

..............

........

. ..... .....

.... . .... .

.. ... ..... ... ...

... ....

:i........ ...

... ...... ... .

.... .. ..

. .....

................i

............. ..

..............

.............

...........

............ .

.............

........... ..

.... ..

. .. . .

.. .........

5.3. Performance Evaluation

5.3.1. Quantification of Errors

In our performance evaluation studies, we studied both registration errors and

volumetric errors. Registration errors are determined by evaluating the Euclidean

distance between the landmark points on the standard image and on the warped

individual image. Two landmarks are used for each structure. One of the landmarks is

chosen on the structure which was aligned, while the other is chosen on a nearby

structure to test the success of alignment in the surrounding area. The anatomic definition

of the four landmarks are as follows:

Cac: Intersection of the cingulate sulcus with the vertical line passing through

anterior commissure. The sagittal slice is chosen to be 1 cm. lateral to the

interhemispheric fissure.

Cs: The superiormost point of the central sulcus. The axial slice is chosen from

the midline sagittal slice where a vertical line drawn on the anterior commissure meets

the hemispheric margin.

Sph: The posteriormost point of the horizontal ramus of the Sylvian fissure. The

sagittal slice is chosen to be midway between the posteriormost point of the insula and

the hemispheric margin, on the axial slice at the level of the anterior commissure.

Spa: The posteriormost point of the ascending ramus of the Sylvian fissure. The

sagittal slice is chosen as in Sph.

Further details on the anatomic definitions can be found in Grachev et al. (1999).

The volume error, Ev is defined to evaluate the goodness of matching even in severe

situations where the structures to be aligned does not even overlap after warping. This is

accomplished by the following calculation.

f max f (x) min fk (x) dx

E =k k (31)

v P(f (x))

iS

where i is an index over all slices in which the structure occurs, j and k are indices over

all subjects, and s is the index of the standard image. The functionfj plots the structure in

slice i, for subject j, and the function P returns the number of points belonging to the

specified structure. This volumetric error definition accentuates the volume difference,

because it also takes into account the type III mismatches. If the features in figure 20 are

studied, the formula in (31) is the same as:

V uV uV

A B Till (32)

V

A

where VTIII denotes the volume of the type III mismatch.In the ideal case, where all

structures in the warped images coincide perfectly, Ev is 1.

In figure 27, the maximum values of the registration and volumetric errors are

shown for the alignment of the cingulate sulcus, and the Sylvian fissure in a set of 9

subjects. The SOMs are chosen to be Id, with 10 nodes for the cingulate sulcus (RCS),

and 2d, with 3x7 nodes for the Sylvian fissure (RSF). As seen in figure 27a, the original

SOM alignment yields more improvement in registration error over the standard linear

Talairach transformation in comparison to regularized SOM alignment. However, we find

the results of RSOM alignment preferable, because it not only improves the registration

error on a scale comparable to SOM alignment, but it also reduces the membrane energy,

by half, as seen in figure 27c. If the volumetric errors are studied, 24-28% improvement

is observed in the alignment of cingulate sulcus, in contrast to 41-59% improvement in

the Sylvian fissure depending on whether regularized SOM or original SOM is used. On

the other hand, the improvement in the registration error is 25-38% and 21-45%

respectively.

It is very important to point out that the error improvements in figure 27 are

contingent upon the selection of the number of SOM nodes. We studied the change of

these errors with respect to the number of nodes on one sample individual image, and

only for the RCS. The results are plotted in figure 28. As seen from figure 28a, as the

number of SOM nodes increases, the quantization error, Eq, decreases. The same trend

applies to the volumetric error and membrane energy. In the volume error graph, K=0 is

used to indicate the initial volume error after the transformation into Talairach

coordinates.

5.3.2. Comparative Evaluation

In this section, the proposed feature based standardization method is compared

with a widely used standardization method, AIR (Woods et al., 1998b). As summarized

in chapter 2, AIR is based on minimizing the global intensity differences between the

given image and the standard image. Minimization is performed through least squared

distance, and warping is achieved through polynomial transformations. For maximum

performance, prior to the application of this algorithm, external tissues, such as skull,

scalp and dura, and the background must be stripped, exposing only the brain. We used

the brain extraction procedure provided in the AFNI package (Cox, 1996) for this

purpose.

O Talairach M Reg. SOM H Orig. SOM

registration error

16

14

12

10-

cs

volume error

14

12Talairach

8

6

4

2

0

Talairach

membrane energy

300000

250000

200000

150000

100000

50000

0

SCing. S. Sylv. F.

Reg. SOM Orig. SOM

0 Original SOM 0 Regularized SOM

Figure 27: Quantification of errors and membrane energy in RSOM and SOM

-7

sOO s02 s03 s04 s06 s07 s08 s09 sl)

4-

25 E

20

15

10

5

K

0

10 15 20 25 30

Volume Error

K

0 5 10 15 20 25 30

35000

30000

25000

20000

15000

10000

5000

Figure 28: Change in quantization error, volume error, and membrane energy

with respect to number of SOM nodes in the alignment ofRCS