Lima J.J.Pedroso, de (ed.). Nuclear Medicine Physics

Подождите немного. Документ загружается.

268 Nuclear Medicine Physics

the large subdivision that is made at the dimension level is to whether to con-

sider only spatial dimensions or to take into account the temporal dimension.

If only spatial dimensions are considered, then it is common to discriminate

the registration methods into 2D/2D, 2D/3D, and 3D/3D.

6.3.3.2 Nature of the Basis of Registration

A distinction is often drawn between approaches that use only a few points

in the images and those that use images as a whole, working on the intensity

of pixels or voxels. We can also distinguish between the extrinsic methods,

which make use of artiﬁcial objects that are placed in the individual(s), and

the intrinsic methods, which are solely based on information contained in the

images. The latter can be further subdivided, depending on what the data are,

into (1) intensity of pixels or voxels, (2) characteristics of images obtained by

segmentation, and (3) characteristic marks that can be seen in the image, for

example, the extremity of the bifurcation of a blood vessel.

6.3.3.3 Nature and Domain of the Transformation

One of the aspects that determines the type of registration algorithm, with

direct implications on its complexity and on the number of variables to

process, is the kind of transformation that is applied and the domain

(global or local) of its application. The specialized literature contains some

splits between algorithms dedicated to rigid transformations and algorithms,

knownas nonrigid, thatintroducedeformations. It is consensual thatnonrigid

methods are used to improve the registration obtained by rigid methods.

In rigid transformations, the correspondence between points in the images

can be described only by translations and rotations, which means that the

relative distance between points is constant. In 2D or 3D registration (the

simplest), only three parameters are needed to characterize the transforma-

tion: two for the translation and one for the rotation. In the registration of two

3D images, three translations and three rotations are necessary.

Rotation and translation are operations that conserve colinearity, which

means that a straight line in one image is transformed into a straight line in

another image. Operations that maintain colinearity and the ratio of distances

(Equation 6.43) between three collinear points, p

, p

, and p

, are called afﬁne

operations.

−p

. (6.43)

Besides rotation and translation, which are special kinds of afﬁne trans-

formations, there are also scaling and skewing transformations. An afﬁne

transformation is often deﬁned as a linear transformation followed by a trans-

lation, that is, an afﬁne transformation transforms a vector x into another

Imaging Methodologies 269

vector x



using the following rule:



= Mx+ t, (6.44)

where M is a combination of rotation, scaling, and skewing and t is a trans-

lation. The elements m

i,j

of the matrix M can have any value. The number of

degrees of freedom increases to 12 in this case, where 6 are due to the rigid-

body transformation, 3 to scaling, and 3 to skewing. Figure 6.33 shows the

effect of applying each of these transformations to a 2D image.

In a rigid transformation, the elements m

of matrix M have to be such that

the determinant is equal to one, which reveals the rotation character of the

matrix. Matrix M is then deﬁned by

M =

⎡

⎣

10 0

0 cos α

−sen α

0 senα

cos α

⎤

⎦

⎡

⎣

cos α

0 −senα

01 0

sen α

0 cos α

⎤

⎦

⎡

⎣

cos α

−sen α

sen α

cos α

001

⎤

⎦

. (6.45)

where the angle α

characterizes the rotation performed around the axis i.

As previously mentioned, a given transformation can be locally or globally

applied, meaning that it can be applied to the whole image or to a part of

the image. The application of local transformations requires special care, as

it may cause violations to local continuity. A single local transformation is,

therefore, rarely applied: Instead, one generally uses several local transfor-

mations in different parts of the image. Figure 6.34 shows some examples of

the transformations described for 2D images.

In a nonrigid transformation, the relative positions of the points are not

invariant, and the number of parameters that characterizes this transforma-

tion is larger than that for an afﬁne transformation. They can, theoretically, be

more than the number of points in the images; in an extreme case, any point

FIGURE 6.33

(a) Original image; (b) change of scale along xx; (c) skewing.

270 Nuclear Medicine Physics

Global

Aﬃne

Nonrigid

Local

FIGURE 6.34

Examples of transformations applied on to 2D images. Local transformations can lead to the

appearance of discontinuities in the image, in case adequate restrictions are not adopted.

(Adapted from P. A. van den Elsen, E.-J. D. Pol, M. A. Viergever. Medical image matching—a

review with classiﬁcation. IEEE Eng. Med. Biol. Mag., 12, 26–39, 1993.)

can move in a particular way. So, these algorithms usually have different

approaches to simplify the problem.

6.3.3.4 Interaction

We can basically distinguish three levels of interaction of the algorithms with

the user.At one end, there are algorithms that do not require any intervention

by the user; and, at the opposite end, there are algorithms that require total

intervention from the user in order to work. The former are known as auto-

matic whereas the latter are known as interactive. Between these extremes,

there are semiautomatic algorithms where the user has an active role in pro-

cessing the registration but does not completely execute it. In these methods,

the user usually participates in the process, whether in the initialization of the

algorithm, in image segmentation, or even by guiding the method through

the acceptance or rejection of hypotheses that are being assumed.

6.3.3.5 Optimization

Optimization involves three different aspects: the measure of similarity that

is used, that is, the objective function that is to be minimized or maximized;

Imaging Methodologies 271

the method for searching for the extremity; and, ﬁnally, the technique for

resampling or interpolation used to translate the image transformation.

6.3.3.5.1 Objective Functions

An exhaustive enumeration of the cost functions that have been used to reg-

ister images is beyond the scope of this book. However, among the functions

we can highlight are the distance between points, difference in intensity,

correlation and mutual information [117].

Markers (intrinsic or extrinsic) are widely used in methods that use

the distance between points. In this case, the problem is mathematically

described by

arg min



−τ(P

)



, (6.46)

where P

and P

represent a set of points to be coregistered in the two

images and τ is the transformation to be determined. Usually



expresses

the Euclidean norm, but other distances can be deﬁned.

In the particular case where the distance between points is used for the

objective function and for a rigid transformation, the algorithm is direct. The

problem, thus, deﬁned is known as “Procrustes” [118] or as “absolute orienta-

tion” in the context of photogrammetry [119–122]. It is known that the solution

is unique if and only if the existence of at least three noncolinear points is guar-

anteed. In 1966, Schönemann [123] proposed a closed-form solution for this

problem that involves two steps [124]: ﬁrst, the translation is determined as a

vector between the centroids of the points in the two images; second, the rota-

tion between the two data sets is calculated (already corrected for translation)

by decomposition of the covariance matrix in singular values. Let ¯μ

and

¯μ

be the centroids of the sets of points considered in the two images. Deﬁn-

ing vectors

and

as the position vectors corrected by the corresponding

centroids, we have

= P

−μ

and

= P

−μ

(6.47)

Thus, the rotation matrix R is such that it veriﬁes the equation

arg min

−R(

)

, (6.48)

which is subjected to the orthogonality constraint

= 1. (6.49)

The rotation matrix is then found by

R = UV

∗

, (6.50)

272 Nuclear Medicine Physics

where matrices U and V* are obtained by decomposition in singular values

of the matrix S =

(

)

−1

S = UΣV

∗

. (6.51)

The methods that are based on the intensity of image values exhibit different

characteristics from those just described. The difference between intensities

[125,126] or, more precisely, the sum of the quadratic pixel-to-pixel differences

(SDQ) is one measure of similarity that can be easily established. Considering

images A and B and a subset Ω of voxels where there is correspondence

between the two images, the sum of quadratic differences can be deﬁned by

SDQ =



X∈Ω

[

A(X) − τ

(

B(X)

)

]

. (6.52)

This is a good objective function for cases where the images differ from

one another only by Gaussian noise [127], but this is a rare situation even

for intramodal registration. Still, the application of this function in magnetic

resonance series is common [128,129].

Correlation is a measure of similarity between two images which assumes

that there is a linear relationship between the intensities of the two coregis-

tered images. The correlation coefﬁcient (CC) is then deﬁned by

CC =



X∈Ω

(A(X) −

A)(τ(B(X)) −





X∈Ω

(A(X) −



X∈Ω

(τ(B(X)) −

. (6.53)

where

A and

B represent the average intensity of images A and B in the same

domain Ω.

Since correlation is not too restrictive relative to the nature of the images, it

can be used between different modalities, as long as there is a linear relation-

ship between voxel intensities, as assumed.Another interesting characteristic

of correlation is the possibility of applying it in the frequency domain, taking

advantage of the use of the Fourier transform [130]. For a rigid transforma-

tion, therefore, a phase difference in the frequency domain will correspond

to a translation in the spatial domain. Rotation can be determined by polar

representation in the frequency domain [131–133].

Mutual information is based on the idea that registration corresponds to the

maximization of the information that is common to the images being regis-

tered. It assumes that the alignment of two images implies the superposition

of the same structures contained in the images and that in the opposite case

there is duplication of the structures in the ﬁnal result. This enables a measure

of information to be used in the registration. The commonest measure for the

information contained in a signal was suggested by Shannon [134] in 1940

and is known as a measure of entropy. The average information, H, given by

Imaging Methodologies 273

a series of symbols j whose probabilities are given by p

, is described by

H =−



log p

. (6.54)

Entropy, H, is maximal if the probabilities p

are all equal, that is, if all the

symbols j occur with the same probability. The entropy is minimal if only one

symbol occurs: In that case, the probability is one for that symbol and zero

for the remaining ones.

Registration implies the existence of two images, so there are two symbols

per voxel (each relative to the intensity of its own image). Consequently,

the concept of joint entropy is deﬁned, which quantiﬁes information in the

two combined images. The more similar the two images, the larger the joint

entropy. The mathematical deﬁnition of the joint entropy between two images

A and B is

H(A, B) =−



(i, j) log p



(i, j). (6.55)

where p

(i, j) represents the joint probability that is usually determined from

the normalized joint histogram. As previously seen, for any superposed voxel

we have a pair of intensities that can be represented on the axes of a chart. The

histogram is then obtained by calculating the relative frequency of each pair.

Figure 6.35 presents some examples for PET images, aligned and misaligned.

The joint histogram is increasingly scattered as the images become more

misaligned, and entropy is used as a measure of this dispersion [135,136].

Mutual information is another measure of the same type used in multimodal

registration and was proposed by both Collignon et al. [137] and Viola and

Wells [127,138]. Pluim [139] reviews three similar deﬁnitions for mutual infor-

mation, of which we only discuss one, as they are all much the same. Mutual

information is directly related to both the entropy of the images and the joint

entropy. Mutual information, I(A, B) of two images A and B is given by

I(A, B) = H(A) + H(B) − H(A, B). (6.56)

FIGURE 6.35

Joint histogram between the same PET image (slice) (a); (b) aligned; (c) translated by 2 pixels; (d)

rotated by 2

◦

274 Nuclear Medicine Physics

The subtractive term points to an equivalence between the maximization of

the mutual information and the minimization of entropy. However, mutual

information presents an advantage over entropy, which is that it also includes

the entropy of the images. This means that, in the optimization process, we

do not obtain a superposition zone containing only the background of the

images, which may occur if only joint entropy is used. Mutual information

can also be thought of as a measure of how well one image explains the other.

6.3.3.5.2 Search for the Minimum

The techniques involved in the search for the minimum of the objective func-

tion directly depend on the characteristics of the function itself. For smooth

and convex functions, Newton’s method and its derivations are usually

preferred, as they are both fast and easy to implement [140–142]. Newton’s

method is an iterative technique to search for zeroes in a function. It relies

on the approximation of a function using the Taylor series expansion [143].

A function, f (x), inﬁnitely differentiable in the neighborhood of a point a,

can be approximated around this point by a polynomial obtained using the

expression

f (x) ≈ f (a) + f



(a)(x −a) +



(a)

(x −a)

+···+

(a)

(x −a)

+···. (6.57)

Considering a linear approximation and assuming that a is a point close to

a zero of function f(x), we deduce the relation

n+1

= x

−

f (x

)



)

, (6.58)

which forms the basis of Newton’s method. From a ﬁrst estimate of a zero

of the function, an increasingly better approximation is obtained. This is a

method that is used for local minima search by searching for the zeroes of the

ﬁrst derivative. Since the objective functions usually have several variables,

Equation 6.58 is adapted to

n+1

= X

−

)

Hf (X

)

−1

∇f (X

), (6.59)

where X

)

...x

is a vector that speciﬁes a point in the space

of function f , ∇ represents the gradient of the function, and H is the Hessian

matrix that contains the second derivatives of the function. Newton’s method

rapidly converges to the solution; but the ﬁrst estimate is crucial, as it must

be close to the solution. Another essential aspect for this method’s success,

mentioned before, is the type of function; for smooth and convex functions,

the method is excellent.

Imaging Methodologies 275

One of the easiest ways of modifying Newton’s method is to include a

variable, γ > 0, to set the step

n+1

= X

−γ

)

Hf (X

)

−1

∇f (X

). (6.60)

The inversion of the Hessian matrix is frequently complex or even impos-

sible. Alternative forms are used in these circumstances, such as the gradient

method (or maximum descent) or conjugated gradient, which do not use

the second derivatives. In the gradient method, the search direction of the

minimum is given by the gradient, so a new point is obtained according to

n+1

= X

−γ∇f (X

). (6.61)

In the case of conjugated gradients, the descending direction is obtained

to speed up the search, using for that purpose a set of directions that are

conjugated between themselves and closest to the gradient [144].

For least squares problems, the Gauss–Newton (GN) algorithm or

Levenberg–Marquardt algorithm [145–147] is also often used, and it com-

bines the advantage of the GN method with that of the gradient. The Gaussian

method can be rapidly derived from Newton’s method, considering its appli-

cation to the mean squares. In the case of the mean squares, the gradient is

given by

∇f (X) = 2J

(X)

f (X). (6.62)

where J

represents the Jacobian of the function. Consequently, the Hessian

is approximately equal to

Hf (p) = 2J

(X)

(X). (6.63)

So, ﬁnally, we obtain

n+1

= X

−



)



−1

)

f (X

). (6.64)

Powell’s method is a frequently used optimization technique [148–151]. It

is characterized by considering a set of mutually conjugated directions and

by searching the minimum in one direction at a time. The appeal of this

method is the absence of the use of gradients. The minimum along a line

is usually determined by means of Brent’s method [152], which combines

a slower but more robust search with a faster search when the function is

smoother. When the function is more difﬁcult, the search for the extreme

involves bisecting the interval. In this case, a triplet of points is considered a,

b, and c so that f(b)<f (a) and f (b)<f (c); and if the function is continuous in

the interval (a, c), then there is a minimum within the interval. This minimum

can, therefore, be isolated by the bisection of intervals (b, c) or (a, b). In fact,

276 Nuclear Medicine Physics

there is no true bisection: given the previous triplet (a, b, c), the new point,

x, tested will be at a fraction of 0.38197 of the largest segment measured

from the central point. This procedure guarantees a golden ratio

∗

between

the points of the triplet, even if it was initiated without its existence; and so,

it assures a linear convergence in the search for the minimum. In situations

where the function is smoother, a faster technique to search for the minimum

can be used, based on a parabolic interpolation. In this case, for the triplet of

abscissae (a, b, c) and ordinates ( f (a), f (b), f (c)), the parabola that contains

them is calculated; and from this, the value of x that is the minimum of the

parabola’s abscissa is found. For noncolinear points we, thus, have

x = b −

(b −a)

)

f (b) −f (c)

−(b −c)

)

f (b) −f (a)

(b −a)

)

f (b) −f (c)

−(b −c)

)

f (b) −f (a)

. (6.65)

After ﬁnding the point x, which is the minimum in a given direction, the

algorithm tests a new direction searching for a new minimum, and so on.

Other methods besides those described above are commonly used in image

registration, in particular genetic algorithms and simulated annealing.

Genetic algorithms are implemented as a computational simulation that

emulates biological evolution to solve a given problem [153–155]. It is, thus,

not uncommon to ﬁnd characteristics such as heredity, mutation, selection,

andrecombinationin genetic algorithms. For a given problem,a setof possible

solutions (population) is given to the algorithm, along with an objective func-

tion that is used to evaluate each particular solution (candidate). The initial

solutions are usually randomly obtained, and it is hoped that the algorithm

will improve them with time. The candidates differ by a set of abstract char-

acteristics that have meaning for the problem being solved. The candidates

are scored and depending on the score are selected to be reproduced, thereby

generating a new population. The candidates of the new population will be

different from the previous ones, despite heredity, due to recombination and

eventual mutation of their characteristics. The stopping criterion of the algo-

rithm can be the number of generations produced or the level established for

the objective function.

Simulated annealing is inspired, as the name implies, from the tempering

of a metal by annealing. This technique involves heating a metal to its melting

point and then cooling it slowly. The heat causes a change in the positions

of atoms, which are in a state of minimum energy, promoting random rear-

rangements of higher energy. The slow cooling facilitates the establishment of

settings with less energy than the original one. The simulated annealing algo-

rithm was independently proposedby Kirkpatrick et al. [156] and Cerny [157],

using the Metropolis algorithm [158]. Similar to the metallurgical process, the

simulated annealing algorithm replaces the current solution by another that

∗

A golden ratio is algebraically deﬁned as a triplet of points

(

a, b, c

)

on a straight line, similar to

those in the relation



ab = ab



bc.

Imaging Methodologies 277

is chosen based on a distribution of probabilities which is a function of a

parameter of temperature. This solution is then accepted or rejected by the

Metropolis condition, which takes into account not only the improvement

of the values of the objective function but also the level of temperature. The

process is repeated until a balance is reached, followed by a new lowering

of temperature and repetition of the process. The fundamental idea of the

algorithm is that at the beginning (high temperature) almost all the states

(solutions) proposed are accepted, including some that worsen the value of

objective function, in order to escape from local minima. As the tempera-

ture falls, the number of possible transitions also decreases until, at the points

when the temperature tends to zero, only the transitions that reduce the value

of the objective function are accepted. At this point, the solution obtained is

the minimum of the objective function.

6.3.3.5.3 Resampling and Interpolation

Since the registration methods are iterative, at each iteration it is necessary to

“move” one of the images, which involves performing a resampling. That is,

it is essential to determine new values of the intensity of pixels or voxels. This

calculation uses interpolation, an approach that, according to Thèvenaz et al.

[159], is the method of recovering a continuous model from discrete values

within a known abscissae range. Therefore, the value, f (X), at any point X of

the space of dimension q is given by linear interpolation, such that

f (X) =



K∈Z

int

(X − K) ∀

X=[x

,...,x

]

∈ M

. (6.66)

Thevalue so determined is, therefore, alinear combination ofknown values,

, sampled in positions K =[k

, k

, ..., k

] and weighted by the values given

by the function ϕ

int

(X − K). The interpolation is, therefore, usually carried out

by convolution with an appropriate kernel. The sync function is an optimal

interpolator; however, in practice, it is difﬁcult to use, because it is spatially

unlimited. Other interpolators are, therefore, preferred, including nearest

neighbor, bilinear, bicubic, and cubic B-splines.

Over time, several methods of interpolation have been proposed that differ

fundamentally and, in practice, in their accuracy and computational complex-

ity [160]. The method to be used must be carefully chosen for each speciﬁc

situation, but bilinear interpolation is one of the methods that gives the best

compromise between accuracy and complexity [115].

6.3.3.6 Modalities

Two kinds of tasks are usually distinguished in relation to modality. Regis-

tration can be done on the same modality, in which case it is monomodal, or

on different modalities, in which case it is multimodal. A typical example of

monomodal registration is what occurs in myocardial SPECT images in which