BSpline registration: Detecting folding

Detect folding

The situation of folding in a 2D or 3D B-spline transformation occurs when there is a negative value in the computed deformation. The computation of the deformation is described in Section [BSplineRegistrationOV_withMath.html#wp1001643 "Deformation" ]. This section describes a method for detecting if such folding occurs given the current positions for the lattice of the B-spline control points. Given the lattice structure of the B-spline control points, folding can be detected by geometric means. In two dimensions, consider any interior control point and its 8 neighboring control points forming a polygon. Construct 8 triangles involving the chosen interior control point each time with the 8 pairings of adjacent control points on the polygon perimeter. The intersection of any one triangle with the other 7 triangles should each be empty in order for there to be no folding. Alternatively, the area of the polygon should be identical to the union of the areas of the 8 triangles.

The issues are similar in three dimensions by considering the 26 neighboring control points forming a polyhedron. In this case, 48 tetrahedrons are constructed involving the chosen interior control point each time along with the triples of control points from the 48 triangles forming the polyhedron mesh. The intersection of any one tetrahedron with the other 47 tetrahedrons should each be empty in order for there to be no folding. Alternatively, the volume of the polyhedron should be identical to the union of the volumes of the 48 tetrahedrons.

A sufficient but not necessary condition to conclude that folding occurs in the case of two dimensions is the premise that any interior control point is outside the polygon formed by its 8 neighboring control points. In the case of three dimensions, the premise is that any interior control point is outside the polyhedron formed by its 26 neighboring control points.

A condition that is both necessary and sufficient to conclude that folding occurs does exist for both 2D and 3D. In the case of two dimensions, again consider the 8 triangles formed by the chosen interior control point each time with the 8 pairings of adjacent control points on the polygon perimeter. It is important that these triangles formed by control points from the lattice be consistently counterclockwise ordered. For each triangle, let A =(a0,a1) be the 2D coordinates of the chosen interior control point, and let B =(b0,b1) and C =(c0,c1) be the 2D coordinates of the other two triangle control points taken in counterclockwise order. The normal to the plane which contains this triangle can be computed from the (right handed) cross product of the two vectors AB and AC. One way to express the calculation of the cross product involves the construction of the following matrix which contains the three triangle point coordinates:

Equation 12

The determinant of this matrix is twice the signed area of the triangle. Starting with consistently counter clockwise ordered triangles from the initially placed lattice of control points, the determinant of this matrix should remain positive as control point positions are moved. The determinant will become negative when the triangle folds.

The premise for the necessary and sufficient condition to detect when folding occurs extends to three dimensions. Consider the 48 tetrahedrons formed by the chosen interior control point each time with the triples of adjacent control points from the 48 triangles forming the polyhedron mesh. It is important that the tetrahedrons formed by the control points from the lattice be consistently counterclockwise ordered, in the 3D sense. For each tetrahedron, let A = (a0,a1,a2) be the 3D coordinates of the chosen interior control point, and let B = (b0,b1,b2), C = (c0,c1,c2), and D = (d0,d1,d2) be the 3D coordinates of the other three tetrahedron control points taken in counterclockwise order. Extending the matrix determinant formulation used for two dimensions to three dimensions, the following matrix is constructed which contains the four tetrahedron point coordinates:

The determinant of this matrix is 6 times the signed volume of the tetrahedron. Starting with consistently counterclockwise ordered tetrahedrons from the initially placed lattice of control points, the determinant of this matrix should remain positive as control point positions are moved. The determinant will become negative when the tetrahedron folds.

Support for color images

The registration of 2D and 3D color images involves first converting the three color channels of values at each sample in the source and target color images to create single channel source and target intensity images. Once the conversion has been done, the single channel intensity source and target images are registered as described in Background. The resulting transformation map is extracted from the registration which is then used to interpolate the original source color image in order to generate the output registered source color image.

Normally, the three color channels represent red, green, and blue intensities, and so the conversion to a single channel value usually represents overall intensity. In this manner, the conversion is simply a weighted average of the three color channels where separate weights can be provided for each channel. The current implementation uses equal weights for each channel.

Support for 2.5D images

A 2.5D image is actually a 3D image containing a series 2D slice images. The registration of a 2.5D intensity or color image involves a 2D registration of each slice in the 2.5D image to a reference slice also in the same 2.5D image. The reference slice can either be a xed slice or the previous slice. The resulting output registered 2.5D image contains the same number of slices as the original 2.5D image where the sample values in a slice of the registered output are from the same slice of the original only "warped" by the B-spline registration transformation.

If the reference slice is fixed, then the reference slice will appear in the registered output as being unchanged. If the reference slice is the previous slice, then the first slice will be registered to the last slice.

Examples of B-spline registration

Figure 1 below shows an example of B-spline 3D image registration. Figure 2 shows the registration deformation image and its lookup table. The registration was done using the default parameter values, which are explained in Section Image types.

Note: that the input and the target image are both images of the same brain, that's why the registration demonstrates a very good overall fitting.

Figure 1. The input source image (A), the input target image (B) and the output registered image (C). Note that images (A) and (B) are both the images of the same brain. That's why it shows a very good overall fitting. Results may vary, if you take different brain images.

Figure 2. The deformation image and its lookup table. The dark areas on the deformation image are related to the areas on images (A) and (B) which did not perform a deformation during the registration or the deformation was relatively small. The light areas are related to the areas on images (A) and (B) which perform a bigger deformation.

Image types

The algorithm can be applied to 2D and 3D color (RGB) and grayscale images.

User Dialogs in MIPAV

Whenever B-spline registration is selected by the user, a dialog is displayed for the user to select the registration reference and the options related to how the registration is to be performed. A different dialog is displayed depending on whether the user selects 2D/3D registration between two images, or 2.5D registration of slices within a single image. The displayed dialogs differ only in how the registration reference is specified which depends on whether 2D/3D registration or 2.5D registration was selected; otherwise, options related to how the registration is to be performed are the same.

The choice of 2D/3D or 2.5D B-spline registration is made by selecting the appropriate menu item from the Algorithms: Registration submenu of MIPAV. Refer to Figure 3 for details.

Figure 3. The choice of 2D/3D or 2.5D B-spline registration is made by selecting the appropriate menu item from the Algorithms: Registration submenu

The 2D/3D B-spline registration is available through the B-spline automatic registration 2D/3D menu item. This menu item is enabled when a 2D or 3D image (intensity or color) is the currently active image. The dialog is displayed only if there is at least one other loaded image in MIPAV which is compatible for registering to the currently active image. A compatible image is one that has the same number of channels of data and the same number of dimensions, but not necessarily the same sizes for each dimension. Refer to Figures 4 and 5. The 2.5D B-spline registration is available through the B-spline automatic registration 2.5D menu item. This menu item is enabled only when a 3D image (intensity or color) is the currently active image. Refer to Figure 6. Since the options related to how the registration is to be performed are identical for 2D/2.5D/3D and intensity/color images, and the only difference is how the registration reference is selected for 2D/3D vs. 2.5D, the dialog for presenting these options to the user is implemented in the same class. This class is JDialogRegistrationBSpline and is in the mipav/model/view/dialogs package. The class dynamically creates a dialog to present to the user with options that are appropriate for the type of registration.

B-Spline Single - Pass Automatic Registration 2D/3D dialog box

Figure 4. B-Spline Single - Pass Automatic Registration 2D/3D dialog box

General
Register	Displays a list of currently active 2D/3D images to registration. The drop down list of options show the names of available images currently loaded to which this image can be registered. This list will only contain the names of images currently loaded in MIPAV which are compatible.	<
Cost function	This drop down list contains the following choices for cost measure to use during registration: Least Squares, Correlation Ratio, and Normalized Mutual Information.
Perform two-pass registra-tion	This dialog is shown with the default options for single-pass registration. If the Perform two-pass registration option is enabled, then the dialog will change to make options for two-pass registration available.
Transformation Options The following is a discussion of the dialog options which control each registration pass, one set of options per pass. For more information refer to Transformation options
Sub-sample image for speed	If this option is enabled, the image will be subsample in power of 2 for each dimension. Subsampling increases performance of the algorithm, but it exhibits lower accuracy of registration.
B-Spline Degree (same for all axes)	This is the degree of the B-Spline to use for all axes of the image. The drop list of options show the values of 1 (linear), 2 (quadratic), 3 (cubic), and 4 which are available. Although, the underlying BSpineBasisf class supports any degree, this dialog restricts the degrees to these values. Choosing a smaller degree executes faster than a larger degree, but a higher degree offers a higher order of continuity.
B-Spline Control Points(same for all axes)	This is the number of B-spline control points to use for all axes of the image.The minimum number of control points is the larger of the number 2 and the degree of the B-Spline. This dialog limits the maximum number of control points to one-half of the number of samples in the smallest dimension.
Gradient Descent Minimize Step Size (sample units)<a name="wp1022724"></a>	This is the size of steps each control point is moved along the gradient-based direction when searching for a minimum of the error. The registration does not take into account the size of an image sample in each dimension, so this step size is specified in units of a sample. E.g., a value of 0.5 will move the control point in increments of half of a sample while searching for the minimum error. The minimum step size is 0.1 sample and the maximum is 5 samples.
Gradient Descent Minimize Max Search Steps	This is the maximum number of steps the control point is moved along the gradient-based direction when searching for a minimum of the error. If the distance from the control point current position to the bounding polygon/polyhedron of its neighboring control points along the gradient-based direction is less than this specified value, then the computed distance value will be used for this control point instead. The minimum number of steps is 1 and the maximum is 100.<
Iteration options (one iteration=move each control point once)
Convergence limit (min change rate for one iteration)		This is a rate of change threshold such that if the rate the error changes between starting and ending of each iteration is smaller than this value, then the registration terminates.
Maximum Number of Iterations		If the convergence limit is not satisfied after this many iterations have been executed, then the registration terminates. The range of possible values is 1 to 100.
Results
Create deforma-tion field image	If this checkbox is enabled, then the deformation image is computed after the last pass of registration is performed. The image that is created has the same dimensions as the target and the registered source images, and has the name of the source image appended with deformation suffix. The deformation image will be automatically displayed along with the registered source image upon completion of registration. Refer to Figure 5 for details.
<b>OK</b>	Applies the algorithm to the input image that you chose in this dialog box.
<b>Cancel</b>	Disregards any changes that you made in this dialog box and closes the dialog box.
<b>Help</b><	Displays online help for this dialog box.</a>

B-Spline Two - Pass Automatic Registration 2D/3D dialog box

If the Perform two-pass registration check box is selected, then the dialog will change to the following:

Figure 5. B-Spline Two - Pass Automatic Registration 2D/3D dialog box

General
Register	Displays a list of currently active 2D/3D images to registration. The drop down list of options show the names of available images currently loaded to which this image can be registered. This list will only contain the names of images currently loaded in MIPAV which are compatible.
Cost function	This drop down list contains the following choices for cost measure to use during registration: Least Squares, Correlation Ratio, and Normalized Mutual Information.
Perform two-pass registra-tion	This dialog is shown with the default options for two-pass registration. If the Perform two-pass registration option is disabled, then the dialog will change to make options for single-pass registration available as shown in Figure 4
Transformation Options The following is a discussion of the dialog options which control each registration pass, one set of options per pass. For more information refer to Transformation options
Sub-sample image for speed	If this option is enabled, the image will be subsample in power of 2 for each dimension. Subsampling increases performance of the algorithm, but it results lower accuracy of registration.
B-Spline Degree (same for all axes)	This is the degree of the B-Spline to use for all axes of the image. The drop list of options show the values of 1 (linear), 2 (quadratic), 3 (cubic), and 4 which are available. Although, the underlying BSpineBasisf class supports any degree, this dialog restricts the degrees to these values. Choosing a smaller degree executes faster than a larger degree, but a higher degree offers a higher order of continuity.
B-Spline Control Points(same for all axes)	This is the number of B-spline control points to use for all axes of the image.The minimum number of control points is the larger of the number 2 and the degree of the B-Spline. This dialog limits the maximum number of control points to one-half of the number of samples in the smallest dimension.
Gradient Descent Minimize Step Size (sample units)	This is the size of steps each control point is moved along the gradient-based direction when searching for a minimum of the error. The registration does not take into account the size of an image sample in each dimension, so this step size is specified in units of a sample. E.g., a value of 0.5 will move the control point in increments of half of a sample while searching for the minimum error. The minimum step size is 0.1 sample and the maximum is 5 samples.
Gradient Descent Minimize Max Search Steps	This is the maximum number of steps the control point is moved along the gradient-based direction when searching for a minimum of the error. If the distance from the control point current position to the bounding polygon/polyhedron of its neighboring control points along the gradient-based direction is less than this specified value, then the computed distance value will be used for this control point instead. The minimum number of steps is 1 and the maximum is 100.
Iteration options (one iteration=move each control point once)
Convergence limit (min change rate for one iteration)		This is a rate of change threshold such that if the rate the error changes between starting and ending of each iteration is smaller than this value, then the registration terminates.
Maximum Number of Iterations		If the convergence limit is not satisfied after this many iterations have been executed, then the registration terminates. The range of possible values is 1 to 100.
Results
Create deforma-tion field image	If this checkbox is enabled, then the deformation image is computed after the last pass of registration is performed. The image that is created has the same dimensions as the target and the registered source images, and has the name of the source image appended with deformation suffix. The deformation image will be automatically displayed along with the registered source image upon completion of registration. Refer to Figure 5for details.
OK	Applies the algorithm to the input image that you chose in this dialog box.
Cancel	Disregards any changes that you made in this dialog box and closes the dialog box.
Help	Displays online help for this dialog box.

This dialog is shown with the default options for a two-pass registration. As previously mentioned and as indicated in the dialog, a single iteration involves the movement of each interior control point in a gradient descent minimization manner.

<a name="wp1003337"></a></p>

B-Spline Single - Pass Automatic Registration 2.5D dialog box

The following is an example of the dialog which is displayed for 2.5D intensity/color image registration: This dialog is show with the default options for a single-pass registration. If the Perform two-pass registration check box is selected, then the dialog will change in a manner similar to that for the 2D/3D image registration shown in Figure 5.

Figure 6. B-Spline Single - Pass Automatic Registration 2.5D dialog box

General
Register	The Register source field is already filled with the name of the currently selected 2.5D image. The reference slice is selected from among the slices in the selected image. Radio buttons are used to select whether the reference slice is the the adjacent (previous) slice or a particular slice number in the image. If the Reference Slice radio button is selected, then a drop down list of slices numbers is enabled where the middle slice is selected by default.
Cost function	This drop down list contains the following choices for cost measure to use during registration: Least Squares, Correlation Ratio, and Normalized Mutual Information.
Perform two-pass registration	This dialog is shown with the default options for single-pass registration. If the Perform two-pass registration option is enabled, then the dialog will change to make options for two-pass registration</em> available. Refer to Figure 7
Transformation Options The following is a discussion of the dialog options which control each registration pass, one set of options per pass. For more information refer to Transformation options.
Sub-sample image for speed	If this option is enabled, the image will be subsample in power of 2 for each dimension. Subsampling increases performance of the algorithm, but it results lower accuracy of registration.
B-Spline Degree (same for all axes)	This is the degree of the B-Spline to use for all axes of the image. The drop list of options show the values of 1 (linear), 2 (quadratic), 3 (cubic), and 4 which are available. Although, the underlying BSpineBasisf class supports any degree, this dialog restricts the degrees to these values. Choosing a smaller degree executes faster than a larger degree, but a higher degree offers a higher order of continuity.
B-Spline Control Points(same for all axes)	This is the number of B-spline control points to use for all axes of the image.The minimum number of control points is the larger of the number 2 and the degree of the B-Spline. This dialog limits the maximum number of control points to one-half of the number of samples in the smallest dimension.
Gradient Descent Minimize Step Size (sample units)	This is the size of steps each control point is moved along the gradient-based direction when searching for a minimum of the error. The registration does not take into account the size of an image sample in each dimension, so this step size is specified in units of a sample. E.g., a value of 0.5 will move the control point in increments of half of a sample while searching for the minimum error. The minimum step size is 0.1 sample and the maximum is 5 samples.
Gradient Descent Minimize Max Search Steps	This is the maximum number of steps the control point is moved along the gradient-based direction when searching for a minimum of the error. If the distance from the control point current position to the bounding polygon/polyhedron of its neighboring control points along the gradient-based direction is less than this specified value, then the computed distance value will be used for this control point instead. The minimum number of steps is 1 and the maximum is 100.
Iteration options (one iteration=move each control point once)
Convergence limit (min change rate for one iteration)		This is a rate of change threshold such that if the rate the error changes between starting and ending of each iteration is smaller than this value, then the registration terminates.
Maximum Number of Iterations		If the convergence limit is not satisfied after this many iterations have been executed, then the registration terminates. The range of possible values is 1 to 100.
Results
Create deformation field image	If this checkbox is enabled, then the deformation image is computed after the last pass of registration is performed. The image that is created has the same dimensions as the target and the registered source images, and has the name of the source image appended with deformation suffix. The deformation image will be automatically displayed along with the registered source image upon completion of registration. Refer to Figure 2 for details.
OK	Applies the algorithm to the input image that you chose in this dialog box.
Cancel	Disregards any changes that you made in this dialog box and closes the dialog box.
Help	Displays online help for this dialog box.

B-Spline Two - Pass Automatic Registration 2.5D dialog box

Figure 7. B-Spline Two - Pass Automatic Registration 2.5D dialog box

General
Register	The Register source field is already filled with the name of the currently selected 2.5D image. The reference slice is selected from among the slices in the selected image. Radio buttons are used to select whether the reference slice is the the adjacent (previous) slice or a particular slice number in the image. If the Reference Slice radio button is selected, then a drop down list of slices numbers is enabled where the middle slice is selected by default.
Cost function	This drop down list contains the following choices for cost measure to use during registration: Least Squares, Correlation Ratio, and Normalized Mutual Information.
Perform two-pass registration<		This dialog is shown with the default options for two-pass registration.
Transformation Options The following is a discussion of the dialog options which control each registration pass, one set of options per pass. For more information refer to Transformation options.
Sub-sample image for speed	If this option is enabled, the image will be subsample in power of 2 for each dimension. Subsampling increases performance of the algorithm, but it results lower accuracy of registration.
B-Spline Degree (same for all axes)	This is the degree of the B-Spline to use for all axes of the image. The drop list of options show the values of 1 (linear), 2 (quadratic), 3 (cubic), and 4 which are available. Although, the underlying BSpineBasisf class supports any degree, this dialog restricts the degrees to these values. Choosing a smaller degree executes faster than a larger degree, but a higher degree offers a higher order of continuity.
B-Spline Control Points(same for all axes)	This is the number of B-spline control points to use for all axes of the image.The minimum number of control points is the larger of the number 2 and the degree of the B-Spline. This dialog limits the maximum number of control points to one-half of the number of samples in the smallest dimension.
Gradient Descent Minimize Step Size (sample units)	This is the size of steps each control point is moved along the gradient-based direction when searching for a minimum of the error. The registration does not take into account the size of an image sample in each dimension, so this step size is specified in units of a sample. E.g., a value of 0.5 will move the control point in increments of half of a sample while searching for the minimum error. The minimum step size is 0.1 sample and the maximum is 5 samples.
Gradient Descent Minimize Max Search Steps	This is the maximum number of steps the control point is moved along the gradient-based direction when searching for a minimum of the error. If the distance from the control point current position to the bounding polygon/polyhedron of its neighboring control points along the gradient-based direction is less than this specified value, then the computed distance value will be used for this control point instead. The minimum number of steps is 1 and the maximum is 100.
Iteration options (one iteration=move each control point once)
Convergence limit (min change rate for one iteration)		This is a rate of change threshold such that if the rate the error changes between starting and ending of each iteration is smaller than this value, then the registration terminates.
Maximum Number of Iterations		If the convergence limit is not satisfied after this many iterations have been executed, then the registration terminates. The range of possible values is 1 to 100.
Results
Create deformation field image	If this checkbox is enabled, then the deformation image is computed after the last pass of registration is performed. The image that is created has the same dimensions as the target and the registered source images, and has the name of the source image appended with deformation suffix. The deformation image will be automatically displayed along with the registered source image upon completion of registration. Refer to Figure 2 for details.
OK	Applies the algorithm to the input image that you chose in this dialog box.
Cancel	Disregards any changes that you made in this dialog box and closes the dialog box.
Help	Displays online help for this dialog box.

Transformation options

The following is a discussion of the dialog options which control each registration pass, one set of options per pass:

B-spline Degree (same for all axes): This is the degree of the B-spline to use for all axes of the image. Even though the BSplineRegistration2Df and BSplineRegistration3Df classes allow for specification of separate degrees for each axis, this dialog restricts the degrees to be the same. The drop list of options show the values of 1 (linear), 2 (quadratic), 3 (cubic), and 4 which are available. Although, the underlying BSpineBasisf class supports any degree, this dialog restricts the degrees to these values.Choosing a smaller degree executes faster than a larger degree, but a higher degree offers a higher order of continuity.

B-spline Control Points (same for all axes): This is the number of B-spline control points to use for all axes of the image. Even though the BSplineRegistration2Df and BSplineRegistration3Df classes allow for specification of separate number of control points for each axis, this dialog restricts the numbers of control points to be the same. Initially, the control points are evenly spaced in each dimension, and the registration algorithm moves the interior control points to minimize the error difference between the target image and the registered source image. The minimum number of control points is the larger of the number 2 and the degree of the B-spline. This dialog option limits the maximum number of control points to one-half of the number of samples in the smallest dimension. Choosing more control points generally results in better overall fitting, but a limitation may be that the control points will be closer to each other and that will limit how much a control point can move to minimize the error (in order to satisfy the constraint that an interior control point must be within the bounding polygon/polyhedron formed by its neighboring control points). Also, since choosing more control points means that control points will be closer to each other, it may be helpful to reduce the maximum number of steps for searching for the gradient descent minimum.

Choosing more control points rather than fewer should result in a longer execution time. Even though more control points means fewer samples (for the same size input image) are affected by a single control point, the fact that there are more control points to move and each one requires a gradient descent sample search for a minimum error most likely results in the longer execution time.

Choosing more control points may require a different convergence limit. Since the convergence limit is tested once per iteration and an iteration involves moving each interior control point once, it may be that more control points could result in more (or less) overall changes in the error between starting and ending one iteration.

Gradient Descent Minimize Step Size (sample size units): This is the size of steps each control point is moved along the gradient-based direction when searching for a minimum of the error. The registration does not take into account the size of an image sample in each dimension, so this step size is specified in units of a sample. For instance, a value of 0.5 will move the control point in increments of half a sample while searching for the minimum error. The minimum step size is 0.1 sample and the maximum is 5 samples.

Gradient Descent Minimize Max Search Steps: This is the maximum number of steps the control point is moved along the gradient-based direction when searching for a minimum of the error. If the distance from the control points current position to the bounding polygon/polyhedron of its neighboring control points along the gradient-based direction is less than this specified value, then the computed distance value will be used for this control point instead. The minimum number of steps is 1 and the maximum is 100.

Convergence limit (min change rate for one iteration): This is a rate of change threshold such that if the rate the error changes between starting and ending of each iteration is smaller than this value, then the registration terminates. As mentioned previously, the rate of change is computed as

(previousError-currentError) / previousError

where previousError is the error measured before starting an iteration and currentError is the error measured after completing an iteration. The error being measured here depends on the particular cost function selected.

Maximum Number of Iterations: If the convergence limit is not satisfied after this many iterations have been executed, then the registration terminates. The range of possible values is 1 to 100.

Output files

After the B-spline registration algorithm finishes running, the following files can be found in the image directory:

• the source and target image,
• the *.xmp files with the recorded history for both source and target images,
• the NonLinear Transformation (*.nlt) file that contains the information about the performed B-spline transformation, such as the dimensionality of the transformed image, the degree of the B-spline, the number of control points, and the values of the control points. This file is used for the Transform Nonlinear algorithm, refer to [TransformNonlinear.html#wp999048 "Transform Nonlinear"],
• the result file.

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