Symmetric Similarity Cost Function

The problem with many image registration techniques is that the image similarity function does not uniquely determine the correspondence between two image volumes. In general, similarity cost functions have many local minima due to the complexity of the images being matched and the dimensionality of the transformation. It is these local minima (ambiguities) that cause the estimated transformation from image T to S to be different from the inverse of the estimated transformation from S to T. In general, this becomes more of a problem as the dimensionality of the transformation increases.

To overcome correspondence ambiguities, the transformations from image $T$ to $S$ and from $S$ to $T$ are jointly estimated. This is accomplished by defining a cost function to measure the shape differences between the deformed image $T \circ h$ and image $S$ and the differences between the deformed image $S \circ g$ and image $T$. Ideally, the transformations $h$ and $g$ should be inverses of one another, i.e., $h = g^{-1}$. In this work, the transformations $h$ and $g$ are estimated by minimizing a cost function

$$C_{\text{SIM}}(T \circ h , S) + C_{\text{SIM}}(S \circ g , T) = \... ..._c(x)) -S_c(x)\vert^2 dx + \int_{\Omega_c} \vert S_c(g_c(x)) - T_c(x)\vert^2 dx$$ (1)

where the intensities of $T_c$ and $S_c$ are assumed to be scaled between 0 and 1. To use this similarity function, the images $T$ and $S$ must correspond to the same imaging modality and they may require preprocessing to equalize the intensities of the image. In practice, MRI images require intensity equalization while CT images do not. A simple but effective method for intensity equalizing MRI data is to compute the histograms of the two images, scale the axis of one histogram so that the gray- and white-matter maximums match, and then apply the intensity scaling to the image.

In practice, the images $S$ and $T$ are discrete and the integrals in Eq. 1 are discretized and implemented as summations:

$$C_{\text{SIM}}(T \circ h , S) + C_{\text{SIM}}(S \circ g , T) = \frac 1 {N_1 N_2 N_3} \sum_{n \in \Omega_d} \vert T_c(h_c(\frac n N)) - S_c(\frac n N)\vert^2 + \vert S_c(g_c(\frac n N)) - T_c(\frac n N)\vert^2$$

$$= \frac 1 {N_1 N_2 N_3} \sum_{n \in \Omega_d} \vert T_c(h_d[n]) - S_d[n]\vert^2 + \vert S_c(g_d[n]) - T_d[n]\vert^2$$

$$= \frac 1 {N_1 N_2 N_3} \sum_{n \in \Omega_d} \vert T_d[N h_d[n]] - S_d[n]\vert^2 + \vert S_d[N g_d[n]] - T_d[n]\vert^2$$ (2)

where the notation $N h_d[n]$ is defined as the $3 \times 1$ column vector $[N_1 h_d^{(1)}[n], N_2 h_d^{(2)}[n], N_3 h_d^{(3)}[n]]^T$ and the terms $T_d[N h_d[n]]$ and $S_d[N g_d[n]]$ are evaluated using trilinear interpolation. Equation 2 is a discrete approximation to the integral cost defined in Eq. 1 which becomes exact as $N_1 \mapsto \infty$, $N_2 \mapsto \infty$, and $N_3 \mapsto \infty$. The discretized cost defined in Eq. 2 demonstrates a trade-off between computational time and accuracy. This trade-off can be exploited by using a low-resolution computational grid $\Omega _d$ during the early iterations of the algorithm when accuracy is less important and using a high-resolution grid $\Omega _d$ in the later iterations to get an accurate final result.

Note that this joint estimation approach applies to both linear and non-linear transformations. In general, the squared-error similarity functions in Eq. 1 can be replaced by any suitable similarity function--mutual information [25,26], demons [6], an intensity variance cost function [27], etc.--where the choice is dependent on the particular registration application (see Discussion).

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