Kluwer - Handbook of Biomedical Image Analysis Vol

26	Farag et al.

The patches that get assigned to the same point on the sphere are aggregated according to the shape categories of the surface components.

The concept of “shape spectrum” features is also included in the COSMOS framework. This allows free-form object views to be grouped in terms of the shape categories of the visible surfaces and their surface areas.

For the recognition purpose, COSMOS adapted a feature representation consisting of the moments computed from the shape spectrum of an object view. This eliminated unlikely object view matches from a model database of views. Once a small subset of likely candidate views has been isolated from the database, a detailed matching scheme that exploits the various components of the COSMOS representation is performed to derive a matching score and to establish view surface patch correspondence.

1.5.1.4 The “Spin Image” Surface Registration

Johnson and Hebert [30] presented an approach for recognition of complex objects in cluttered 3D scenes. Their surface representation, the “spin” image, comprises descriptive images associated with oriented points on the surface. Using a single point basis, the positions of the other points on the surface are described by two parameters. These parameters are accumulated for many points on the surface and result in an image at each oriented point which is invariant to rigid transformation.

Through correlation of images, point correspondences between a model and scene data are established. Geometric consistency is used to group the correspondences from which plausible rigid transformations that align the model with the scene are calculated. The transformations are then refined and verified using a modified ICP algorithm.

The spin image is generated by first considering an oriented point (a 3D point with a normal direction) and defining its basis. The basis is defined using the tangent plane perpendicular to the point direction and the point direction itself. A spin-map is then defined using the point basis. In this spin map, any other point on the surface is related to the oriented point by two parameters, one is the perpendicular distance to the oriented point line direction and the other is the signed perpendicular distance to the plane passing through the oriented point and perpendicular to the point direction.

Medical Image Registration

27

Figure 1.8: The spin image generation process can be visualized as a sheet spinning around the oriented point basis, accumulating points as it sweeps space.

The term “spin-map” comes from the cylindrical symmetry of the oriented point basis; the basis can spin about its axis with no effect on the coordinates with respect to the basis. To create a spin image, first an oriented point on the surface is selected. Then for each other point on the surface, the spin-map parameters are computed. These parameters are then accumulated in a 2D array. Once all the points on the surface have been processed, the 2D array is converted into a gray image. Figure 1.8 shows the spin-image generation process and visualizes it as a sheet spinning around the oriented point basis, accumulating points as it sweeps space. Figure 1.9 shows examples of spin images generated for the surface of a statue. The darker the pixel, the higher the number of points projected into this location.

They compared spin images using linear correlation coefficients. They used the magnitude of the correlation coefficients as well as the confidence in the correlation results which is measured by the variance of the correlation coefficient. They also modeled the effect of clutter and occlusion to predict a lower bound on the correlation coefficient.

28							Farag et al.

Figure 1.9: Examples of spin images generated for the surface of a statue. The darker the pixel, the higher the number of points projected into this location.

1.5.1.5 The “Surface Signature” Surface Registration

The surface signature algorithm captures the 3D curvature information of any free-form surface and encodes it into a 2D image corresponding to a certain point on the surface. This image is unique for this point and is independent from the object translation or orientation in space.

The process starts by identifying special points on the surface. These points are called Important points due to the information they carry. Then an image capturing the surface curvature information seen from each important point, is formed. This image is called the Surface Signature at this point because it is almost unique for each point on the surface. Surface registration is then performed by matching signature images of different surfaces and hence finding

Medical Image Registration

29

corresponding points in each surface. A point P on a surface/curve S, is called important point, PI , if and only if the absolute value of the curvature at this point is larger than a certain positive value (a threshold).

A = {PI } = {P S| |Curv(P)| > , > 0}

(1.15)

As the important points are landmarks, one may expect that they are stable for the same object. However, due to scanning artifacts, their number and locations may vary. By adjusting the curvature threshold, a common subset can be found. Otherwise, the object has either suffered from non-rigid transformations or its visible surface has no landmarks.

The signature, computed at each important point, encodes the surface curvature seen from this point, thus giving it discriminating power. As shown in Fig. 1.10 (a), for each important point P A defined by its 3D coordinates and the normal UP at the patch where P is the center of gravity, each other point Pi on the surface can be related to P by two parameters: The distance

		di = ||P − Pi||				(1.16)
	and the angle
	α	cos		1	UP · (P − Pi)	(1.17)
			−		||P − Pi||
	i =

Also you can notice that there is a missing degree of freedom in this representation which is the cylindrical angular parameter. At each location in the image, the gray value encodes the angle

βi = cos−1(UP · UPi )

(1.18)

This represents the change in the normal at the surface point Pi relative to the normal at P. Figure 1.10 (b) shows some signature images taken at different important points on the skin model of a patient’s head.

The next step in the registration process is to match corresponding signature images of two surfaces. The ultimate goal of the matching process is to find at least a three-points correspondence in order to calculate the transformation parameters. The benefit of using the signature images to find the correspondence is that we can now use image processing tools in the matching, hence reducing the time taken to find accurate transformation. One such tool is Template Matching in which a value defines how well a portion of an image matched a template.

30	Farag et al.

αi	α

UP UP

βi

αi

	P	di	βi
UPi

di

Pi

UPi

d

SPS image at P

Figure 1.10: (a) For each important point on the surface a signature image is obtained. This image encodes the angle difference between the normal at the point in focus and the normal at each other point on the surface. (b) Examples of signature images taken at different important points locations. Notice how the image provides features the curvature information. The dark intensity in the image represents a high curvature seen from the point while the light intensity represents a low curvature. Also notice how different the image corresponding to a location is from images of other locations.

Medical Image Registration

31

The end result of the matching process is a list of groups of likely threepoints correspondences that satisfies the geometric consistency constraint. The list is sorted such that correspondences that are far apart are at the top of the list. A rigid transformation is calculated for each group of correspondences and a verification stage [9] is performed to obtain the best group. Detailed discussion concerning the surface signature sensitivity and robustness can be found in [26].

1.5.2Maximization of Mutual Information (MI) Algorithm

MI is a basic concept from information theory, measuring the statistical dependence between two random variables or the amount of information that one variable contains about the other. The MI registration criterion used states that the MI of corresponding voxel pairs is maximal if the two volumes are geometrically aligned [31]. No assumptions are made regarding the nature of the relation between the image intensities in either modality.

Consider the two medical volumes to be registered as the reference volume

R and the floating volume F. A voxel of the reference volume is denoted R(x), where x is the coordinates vector of the voxel. A voxel of the floating volume is denoted similarly as F(x). Given that T is a transformation matrix from the coordinate space of the reference volume to the floating volume, F(T (x)) is the floating volume voxel associated with the reference volume voxel R(x).

MI seeks an estimate of the transformation matrix that registers the reference volume R and floating volume F by maximizing their mutual information. The vector x is treated as a random variable over coordinate locations in the reference volume. Mutual information is defined in terms of entropy in the following way [25]:

I(R(x), F(T (x))) ≡ h(R(x)) + h(F(T (x))) − h(R(x), F(T (x))).

(1.19)

where h(R(x)) and h(F(T (x))) are the entropy of R and F, respectively. h(R(x), F(T (x))) is the joint entropy. Entropy can be interpreted as a measure of uncertainty, variability, or complexity. The mutual information defined in Eq. (1.19) has three components. The first term on the right is the entropy in the reference volume, and is not a function of T . The second term is the entropy of the part of the floating volume into which the reference volume projects. It

p(r(x), F(T (x))) · log2

32	Farag et al.

encourages transformations that project R into complex parts of F. The third term, the (negative) joint entropy of R and F, contributes when R and F are functionally related. Maximizing MI tends to find as much of the complexity that is in the separate volumes (maximizing the first two terms) as possible so that at the same time they explain each other well (minimizing the third term) [25, 31].

1.5.2.1 Computation of MI Metric

Based on the definition of relative entropy, also known as Kullbak Leibler distance, between two probability mass functions, Eq. (1.19) can be written in terms of probability distribution functions as follows

I(R(x), F(T (x)))

≡

R(x),F(T (x))

p(r(x), F(T (x)))

. (1.20)

pR(x)(r(x) · pF(T (x))(F(T (x)))

where p(x, y) is the joint distribution function and px(x) is the marginal probability mass functions. The marginals can be obtained directly from the joint probability function. The joint probability mass function p(r(x), F(T (x))) will be approximated by the normalized joint histogram H(r(x), F(T (x))). Here, normalization refers to scaling of the histogram, such that the sum of approximated probabilities equals 1.0. The marginals are then approximated from H( ) by summation over the rows, and then the columns. Computation of H( ) involves a complete iteration over each sample in the floating volume. For each sample, the transformation T is applied, to arrive at a coordinate set in the image coordinate system of the reference volume. If the transformed coordinate is outside the measured reference volume, then the remaining operations are not executed, and the process starts again with the next sample in the floating volume. Otherwise, a sample in the reference volume at the transformed coordinates is approximated using trilinear interpolation, and discretized. The two samples, one from the floating volume, and one from the reference volume, are then binned in the joint histogram.

Computation of the joint histogram involves the processing of each sample in the floating volume, application of a transformation to the coordinate of the sample in the floating volume to obtain a coordinate in the reference volume, interpolation in the reference volume, and binning in the joint histogram. For a typical 256 by 256 by 20 MRI volume, there are thus 256 by 256 by 20 = 1,310,720

Medical Image Registration

33

samples to process. Following computation of the joint histogram, normalization and computation of the marginal histograms must be performed. This involves one pass over the joint histogram, therefore processing 256 by 256 = 65,536 elements from the joint histogram. This processing consists of normalization, and summation to compute the marginal histograms. Following this operation, the mutual information metric itself may be computed. This processing involves computation of the sum given in Eq. (1.19), and involves one pass over the joint histogram, therefore processing 256 by 256 = 65,536 elements that compose the sum given in Eq. (1.19). Therefore, computation of the joint histogram is by far the most computationally costly component in computation of the mutual information metric. Therefore, performance may be best increased by decreasing the execution time of the computation of the joint histogram. Computation of the joint histogram is an amenable problem for parallel execution, as computation of a part of the joint histogram does not depend on the computational results of any other part of the joint histogram, allowing individual bins, or entire regions of the joint histogram to be computed independently, and then merged to form the total joint histogram. Figure 1.11 illustrates such parallel architecture.

														1
																			256
									256
					1
													2
																		256
																				Σ
									256
				2																		256
																								256
												n					256
										256
					n
volume		sub-volumes						sub-joint histograms															joint histogram

task 1		task 2

Figure 1.11: Illustration of the parallel computation of the joint histogram nec-

essary for computation of the mutual information metric.

34	Farag et al.

1.5.2.2 Optimization Methods

The mutual information metric provides a quantitative measure of spatial alignment between two volumes, given a choice of registration parameters. To obtain the best alignment, it is necessary to maximize the metric. Maximization of the metric, which is parameterized in terms of the registration parameters, is numerically accomplished with the use of an optimization algorithm. In the original formulation of registration by maximization of mutual information in [31], Powell’s multidimensional optimization method, with Brent line minimizations was used for maximization of the mutual information metric [32]. Subsequently, [33] compares different classical optimization methods maximizing the mutual information metric. One such method included was the use of the classic Nelder and Mead or simplex algorithm for maximization. This method solely uses the objective function directly for optimization, and therefore does not require the expensive computation of derivatives. This method is a geometry-based method, using the geometric operations of contraction, expansion, and reflection to manipulate a simplex to a maximum of the objective function. GA was also used as an optimization criteria in maximizing the mutual information metric as demonstrated in [34].

1.6 Practical Examples

1.6.1 Composite Approach for Mutlimodal Registration

A composite registration approach can perform a fast, six degrees of freedom registration and accurately locate the position and orientation of medical volumes (obtained from CT/MRI scans for the same patient) with respect to each other. The technique uses surface registration technique followed by a volume registration approach. The advantage of this combination is to have an accurate alignment and to reduce the time needed for registration. The surface registration uses the surface signature approach and for the volume registration, the maximization of mutual information (MI) is used as a matching criterion and is enhanced by a genetic based optimization technique. Figure 1.12 shows a block diagram of the composite registration. Figure 1.13 shows some results and Table 1.2 shows approximate registration time for both single and composition registration.

Medical Image Registration				35
volume 1
		Initial
	Surface			Maximizing
	Registration	estimate		MI using GA

volume 2

Figure 1.12: Two different volumes are fed into the system. Surface registration provides an estimate transformation between them. This result is further enhanced by maximizing the mutual information. A Genetic Algorithms technique is used in the optimization process.

1.6.2 Multiple Sclerosis (MS) Application

Multiple sclerosis is a disease of the human central nervous system, affecting approximately 250,000 to 350,000 people in the United States alone [35]. MS results

BeforeRegistration

After Registration

Figure 1.13: Two different signature samples. The upper two images are from the MR dataset; while the lower two are from the CT dataset. This figure shows similarity of the signature images of each matched pair, although they are produced from two different datasets (CT and MR) with a different 3D volume size (number of slices).