Consistent Landmark and Intensity-based Registration

The consistent landmark and intensity image registration algorithm is outlined in Fig. 1. It is assumed that the images being registered have been rigidly rotated and translated to put them into a standard orientation, such as the Talairach coordinate system[35], before applying this procedure. In this work the data sets did not need to be scaled to the same size before applying the algorithms since the landmark registration takes care of the global scale differences. The first step of the algorithm is to produce a good initial nonrigid registration using a landmark initialization step. This step consists of (1) picking corresponding landmarks in the two images, (2) solving the unidirectional TPS algorithm modified to produce periodic boundary conditions for the forward and reverse transformations, and (3) averaging the forward transformation with the inverse of the reverse transformation, and vice versa. The full details of the landmark initialization are described in Section II-B.

Figure 1: Flow chart describing the steps of the Consistent Landmark and Intensity Image Registration algorithm. The images being registered are assumed to be rigidly aligned before starting this procedure.

After the landmark initialization step, the iterative consistent landmark thin-plate spline algorithm described in Section II-C is used to produce a consistent set of forward and reverse transformations. This algorithm jointly estimates a set of transformations that minimize both the inverse consistency error and the bending energy of the thin-plate spline model while maintaining exact correspondence at the landmarks. The resulting forward and reverse transformations have orders of magnitude less inverse consistency error than the original unidirectional TPS transformations as shown by the experiments in the Results Section.

The last step of the algorithm is to use the consistent intensity registration algorithm[32,33,34], that is briefly described in section II-D, to refine the transformations based on matching the intensities of the images. This step matches the images in regions away from the landmarks by minimizing the intensity differences in these regions. The intensity matching does little in regions near corresponding landmarks since these regions have similar intensity patterns that have all ready been matched by the landmark registration. During the intensity matching step, the landmark correspondence error increases in regions where there are bad landmark initializations. The landmark correspondence error also increases in this step due to the TPS regularization model that can pull corresponding landmarks apart. The landmark registration error can be minimized by applying the consistent landmark registration step followed again by the intensity registration step.

The process of alternating between matching the landmarks and then the image intensities is repeated until an appropriate stopping criteria is met. In this work, a fixed number of iterations was used as the stopping criteria. Alternatively, the algorithm could be stopped after an acceptable intensity similarity and landmark error cost is achieved. The optimal strategy for stopping the algorithm can be quite complex and will be studied in future work.

The consistent landmark and intensity-based TPS algorithm can be thought of as estimating a consistent set of forward and reverse transformations that minimize the intensity differences between two images while being guided by the landmark correspondences. The landmarks guide the solution by initializing the consistent intensity registration algorithm with transformations that are nearly inverse consistent and have exact correspondence at the landmarks. This initialization helps the consistent intensity registration avoid some local minima and therefore produce more biologically relevant correspondence maps.

The final registration is determined from the intensity information alone. The landmark matching is used to get the two images close in a global sense and then the intensity information is used to fine tune the registration in the neighborhood of the landmarks. The choice not to use the landmarks for the final registration was motivated by the fact that the landmarks are generally located on object edges or in regions of changing intensity in the image. Therefore once the landmark matching gets corresponding regions of the images close, the intensity matching finishes the job by matching all of the points in the neighborhood of the landmark based on their intensities.

Another reason that the final registration only used intensity information was that we encountered problems when we tried to force the landmarks and intensity to match at the same time. The problems occured in regions where the landmarks did not correspond to the same location on the intensity profile in both data sets. The Jacobian of the transformation would go negative at the landmark locations when the intensity and landmark forces moved in different directions. This problem was due to the fact that the landmark driving force was focused at the landmark and the intensity driving force was distributed. The Jacobian of the transformation would go negative as the landmark point moved in one direction and the material surrounding the landmark would go in the opposite direction.

The following notation will be used throughout the rest of the paper. The variables $q_i$ and $p_i$, for $i=1, \dots, M$, denote the $M$ corresponding landmarks in the template $T$ and target $S$ images, respectively. The domain of the template image $T$ and target image $S$ is denoted by $\Omega$. The forward transformation $h:\Omega \rightarrow \Omega$ is defined as the mapping that transforms $T$ into the shape of $S$ and the reverse transformation $g:\Omega \rightarrow \Omega$ is defined as the mapping that transforms $S$ into the shape of $T$. The forward and reverse displacement displacement fields are defined as $u(x)=h(x)-x$ and $w(x)=g(x)-x$, respectively. The inverse of the forward and reverse transformations denoted by $h^{-1}(x)$ and $g^{-1}(x)$, respectively, can be expressed in terms of the displacement fields $\tilde{u}(x)= h^{-1}(x)-x$ and $\tilde{w}(x)=g^{-1}(x)-x$, respectively.

Copyright © 2002 • The University of Iowa. All rights reserved. Iowa City, Iowa 52242
Questions or Comments: gary-christensen@uiowa.edu