Difference between revisions of "B-Spline Automatic Registration"

Latest revision as of 18:14, 9 May 2012

This chapter describes the B-spline based registration of images in three dimensions (3D). The description specializes to two dimensions (2D) is in the obvious manner. Chapter is continued on the following pages:

• B-Spline Automatic Registration - you are here

• BSpline registration: Detecting folding

• Examples of BSpline registration

• B-Spline Automatic Registration Dialog

Class implementation for MIPAV

The following classes are part of the mipav/model/structures package:

• B-splineBasis
• B-splineBasisDiscretef
• BSplineLattice2Df
• BSplineLattice3Df
• BSplineRegistrationBasef
• fBSplineRegistration2Df

The following classes are also part of the mipav/model/algorithms/registration package:

• AlgorithmRegBSpline
• AlgorithmRegBSpline2D
• AlgorithmRegBSpline3D
• AlgorithmRegBSpline25D

Background

Let the target image have dimensions <math> (p_o)\,(p_o)\,(p_o) </math> with image indices satisfying <math> 0\le j_0<p_0 </math>, <math> 0\le j_1<p_1 </math>, and <math> 0\le j_2<p_2 </math>. Let the source image have dimensions <math> (q_o)\,(q_o)\,(q_o) </math> with image indices satisfying <math> 0\le k_0<q_0 </math>, <math> 0\le k_1<q_1 </math>, and <math> 0\le k_2<q_2 </math>. The B-spline volume function maps indices <math> j_0,j_1,j_2 </math> via equation 1.

where the <math> N_{i_e, d_e}(u) </math> functions are the B-spline basis functions, and the B values are the control points. The numbers <math>

n_0,n_1,</math> and  <math> n_2 </math>, are user selectable, the degrees 

<math>

d_0,d_1,</math> and  <math> d_2

</math> for the basis functions are user selectable, and the B-spline basis functions are open and uniform.

The control points have the following characteristics:

• Each control point has coordinates from the unit cube. In two dimensions, each control point has coordinates from the unit square.
• The initial placement of the control points depends on the B-spline basis and must yield the identity function; that is,

<math>

M_{j_0j_1j_2} = (\frac {j_0} {p_0-1}), (\frac{j_1} {p_1-1}), (\frac {J_2)} {P_2+1})

</math>. The initial placement of the control points is described in Section "Improving performance".

• A control point is said to be in the outer layer if one of its indexes is either the minimum index (zero) or the maximum index (one less than the number of control points in a given dimension). The position of the control points in the outer layer are fixed as follows:

This constraint implies that the <math>

0\le k^\prime_w<q_w </math>' 'for all w.

• A control point not in the outer layer is said to be an interior control point. The position of the interior control points is constrained such that its movement must not cause any folding in the B-spline transformation. Folding is defined to occur when there is a negative value anywhere in the computed deformation for the current B-spline transformation. The computation of the deformation is described in Section "Deformation". The method by which folding is detected is described in Section "Detect folding".

Each <math>

M_(j_0)(j_1)(j_2)

</math> contains normalized 3D coordinates (i.e., in the unit cube) that identify a location in the source image corresponding to the sample at <math> (j_0,j_1,j_2) </math> in the target image. The <math>

M_(j_0)(j_1)(j_2)

</math> coordinates are related to sample coordinates in the source image by <math>

M_{j_{0}, j_{1}, j_{2}} =\left( \frac {k^\prime_0}{q_0-1}, \frac{k^\prime_1}{q_1-1}, \frac{k^\prime_2}{q_2-1}\right)

</math> where the output indices <math> k^\prime_0, k^\prime_1, k^\prime_2 </math> are considered to be continuous values; that is, they are not necessarily integer valued. If the source image's samples are <math> s_k\_0,k_1,k_2 </math>, where <math> 0\le k^\prime_w < q_w </math> for all w , we need to evaluate the image at the nonintegral values <math> k^\prime_0, k^\prime_1, k^\prime_2 </math>. We can use trilinear interpolation to do this (use bilinear interpolation in the case of two-dimensional (2D) registration). Let <math> k_w= |k^\prime_w| </math>, the largest integer smaller than <math> k^\prime_w </math>. Define <math> \Delta w = k^\prime_w - k_w </math>. The interpolated image value is:

Equation 3
<math> U_0 = (1-\Delta_0)(S_{k_0,k_1,k_2}) + (\Delta_0)(S_{k_0 + 1,k_1,k_2}) </math>

<math> U_1 = (1-\Delta_0)(S_{k_0,k_1+1,k_2} + (\Delta_0)(S_{k_0,k_1 + 1,k_2}) </math>

<math> U_2 = (1-\Delta_0)(S_{k_0,k_1,k_2+1} + (\Delta_0)(S_{k_0+1,k_1,k_2+1}) </math>

<math> U_3 = (1-\Delta_0)(S_{k_0,k_1+1,k_2+1} + (\Delta_0)(S_{k_0+1,k_1+1,k_2+1}) </math>

<math> \nu_0 = (1 - \Delta_1)U_0 + (\Delta_1)(U_1) </math>

<math> \nu_1 = (1 - \Delta_1)U_2 + (\Delta_1)(U_3) </math>

<math> S(\kappa_0, \kappa_1, \kappa_2) = (1- \Delta_2)\nu_0 + (\Delta_2)(\nu_1) </math>

If the target image samples are <math> T_{j_0, j_1, j_2} </math>, where <math> 0 \le j_w < p_w </math> for all w, then the source image S is trilinearly interpolated at <math> (k^\prime_0,k^\prime_1, k^\prime_2) </math>, which are the coordinates of <math>

M_(j_0)(j_1)(j_2)

</math>. This means that the currently registered source image, denoted , has the same dimensions as the target image and can be defined in terms of the target image samples as <math> S^\prime(j_0,j_1,j_2) = S(k^\prime_0,k^\prime_1, k^\prime_2) </math>. A measure of error between the two images Sj0,j1,j2 and Tj0,j1 ,j2 can be computed. A number of possible error measure functions are defined in Section 1.2. The goal is to select the control point positions to make the error measure as small as possible. This is done by moving a single control point at a time using gradient descent.

The gradient of the error function is estimated for the control point Bi0,i1,i2 by means of finite differences, where changes in the error function are measured for small steps of the control point from its current position in each of the coordinate axis directions. Let , and be small step sizes equivalent to one sample increment in the target image,

The gradient of the error function is estimated at each control point Bi0,i1,i2. The gradient at each control point is denoted by

The gradient descent step is as follows. The current control point under consideration is Bi0,i1,i2. The ray originating at this point and having direction of greatest decrease in E is

where the direction is the unit length vector

and where is computed using equation 5. The idea is to choose t > 0 to minimize

A calculus approach which sets the first derivative to zero and solves for t yields an expression for e(t) that is really complicated. Instead, you should choose a set of t-values, t'1 through t'm, evaluate e(t'i') and find the minimum value. Suppose this occurs at t'k. Replace this in equation 6 to obtain the new location for the control point.

Only those t-values which do not result in the occurrence of folding in the B-spline transformation are considered. The method described in Section "Detect folding" is used to determine if folding occurs for a given set of control point positions. Note that sample t0 (i.e., t=0) is equivalent to the location of the control point before any movement attempt, and by definition, sample t0 does not cause folding. A test is made to determine whether folding occurs for tm, the extreme sample in the set of chosen t-values. If folding does not occur at tm, then folding does not occur at any of the other t-values between t0 and tm. If folding does occur at tm, then a bisection method is employed to find the tf old value between t0 and tm where folding begins to occur. The search for a minimum value of e(t) by sampling must be limited to the t < t'f old.

Error measures

All of the error measures between two images defined in this section assume that the two images have the same dimensions. This is the case for the target image T (j0,j1,j2) and registered source image S (j0,j1,j2) = S(k0,k1,k2).

The particular error measures defined here are only concerned with corresponding samples from the two images or histograms of the values of the samples from the two images, and are not concerned with whether the images have two or three dimensions. For simplicity, each image measure will be defined with a single index to access its samples.

For more information, refer to Cost functions used in MIPAV algorithms.

Least squares measure

The least squares measure is defined by

where the summation is over all the indices for the target image voxels.

Here, S'`(j) is a registered source image and T(j) is a target image. The minimum value of this measure is zero; the maximum value depends on the range of values in the target and registered source images.

For more information, refer to Cost functions used in MIPAV algorithms.

Correlation ratio measure

See also: Cost functions used in MIPAV algorithms. The correlation ratio measure is defined by

where

• Var(S) is the variance of the registered source image S.
• k is the number of intensity histogram bins equally distributed over the range of intensity values of the target image T.
• S'k'' is the k-th isoset defined as the set of intensities in registered source image S' at positions where the intensity in target image T is in the k-th intensity bin.
• nk is the number of elements in the set Sk, such that N =Sk ''n'k.
• N is the number of samples in the target image T which is the same number of samples in the registered source image S'.
• Var(Sk) is the variance of the values in the set Sk.

The range of values for this measure is [0,1].

Normalized mutual information measure

The normalized mutual information measure is defined by

where

is the standard entropy definition for probability pij estimated using the (i,j) joint histogram bins, and similarly for the marginals H(S') and H(T).

The standard entropy calculation

for an image X of N samples and user-selectable k histogram intensity bins equally distributed over the range of samples in image X can be rewritten as

where nk is the histogram value for the k-th intensity bin and N =S k nk. The computational effort can be reduced by factoring to get

where log n'k can be a precomputed lookup table of logarithms of integer values.

The range of values for this measure is [0,1].

For more information, refer to Cost functions used in MIPAV algorithms.

Improving performance

The gradient descent applied to a control point at a time can be quite slow. Performance can be improved by taking advantage of the local control aspect of using B-splines. Since only one control point is being moved at a time, the local control aspect of using B-splines means that only some of the registered source image values will be different. Note that for an image of a given size, the number of image values, under local control, affected by moving a control point decreases as the number of control points increases.

One performance improvement can be achieved by precomputing the B-spline basis function values once for each possible sample for each dimension. Since the registration occurs by transforming the source image into the space of the target image, all samples for B-spline basis function evaluation always occur at the discrete sample coordinates for the target image. Since the number of target images samples for a dimension is xed, the B-spline basis function will be evaluated on the [0, 1] interval at equally spaced intervals (we are using open, uniform B-splines). For each dimension, a 2D array is created which stores B-splines basis function evaluations for each sample with the weight associated with each control point. Because of local control, many of these weights will be zero indicating that a sample is not affected by a particular control point. In the process of creating this 2D array of precomputed basis function evaluations, the index range of samples affected by each control point can be determined. The algorithm to determine this range does not examine the weights as they are computed; instead, the range of control points affecting a single sample are determined from the computed knot index and the known B-spline degree.

Another performance improvement involves maintaining intermediate error measure calculations. For example, in the case of the least squares error measure, an image is maintained of the current error calculations, that is, the target image subtracted from the current registered source image. A cumulative total error value is also maintained, which is the sum of the squared error values for all current values in the error image. When a control point is moved, the old values that are being replaced are first subtracted from the total error, and then the new values to be stored are added to the total error. Similar approaches are performed for the correlation ratio and normalized mutual information error measures such as storing cumulative histogram values and images which contain predetermined histogram bin indexes.

Initial placement of control points for identity map

Referring to equation 1 for the definition of the B-spline mapping function, the identity map is the one which yields

M_{j_0j_1j_2} = (\frac {j_0} {p_0-1}), (\frac{j_1} {p_1-1}), (\frac {J_2)} {P_2+1})