Estimation Procedure

The Fourier series parameterization in Eq. 6 is a multi-resolution decomposition of the displacement fields. Let $\Omega_d[r] = \Omega_d \backslash G[r]$ represent a family of subsets of $\Omega _d$ where $G[r] = \{ n \in \Omega_d \vert r_1 < n_1 < N_1 - r_1; r_2 < n_2 < N_2 - r_2; r_3 < n_3 < N_3 - r_3\}$ and the set subtraction notation $A \backslash B$ is defined as all elements of A not in B. Figure 2 illustrates the definition of $\Omega_d[r]$ using a 2D example. In practice, the low frequency basis coefficients are estimated before the higher ones allowing the global image features to be registered before the local features. This is accomplished by replacing Eq. 7 by

$$u_d [n,r] = \sum_{k \in \Omega_d[r]} \mu[k] e^{j<n,\theta[k]>} \quad w_d [n,r] = \sum_{k \in \Omega_d[r]} \eta[k] e^{j<n,\theta[k]>} .$$ (11)

where $r \in \Omega_d$ determines the number of harmonics used to represent the displacement fields. Let the set of multi-resolution transformations be defined as $h_d[n,r] = \frac n N + u_d[n,r]$ and $g_d[n,r] = \frac n N + w_d[n,r]$. The components of $r$ are initially set small and are periodically increased throughout the iterative minimization procedure. The set $\Omega_d[r]$ can be replaced by $\Omega _d$ when all of the components of $r$ are greater than or equal to $(N -1)/2$ since the set $G[r]$ is empty. The constants $r_1$, $r_2$, and $r_3$ represent the largest $x_1$, $x_2$, and $x_3$ harmonic components of the displacement fields. Each displacement field in Eq. 11 is efficiently computed using three $N_1 \times N_2 \times N_3$ FFTs, i.e., each component of the $3 \times 1$ vectors $u_d$ and $w_d$ are computed with a FFT after zeroing out the coefficients not present in the summations.

Figure 2: Example of elements contained in an $8 \times 8$ lattice. (a) $\Omega _d$ contains all of the points while $G[2,2]$ contains the points within the doted lines; (b) Elements of $\Omega_d[2,2] = \Omega_d \backslash G[2,2]$.

The image registration problem can be stated mathematically by combining Eqs. (2), (3), (4), and (9) and estimating the set of parameters $\{ \hat{\mu}[k],\hat{\eta}[k] \}$ that minimize the combined cost function

$$\begin{multline} C(\mu,\eta,r) = \sigma \frac 1 {N_1 N_2 N_3} \sum_{n \in \Ome... ...\mu^{\dag }[k] D^2 [k] \mu [k] + \eta^{\dag }[k] D^2[k] \eta[k]. \end{multline}$$

The constants $\sigma$, $\chi$ and $\rho$ are used to enforce/balance the constraints. The effects of varying these parameters are studied later in the paper.

The basis coefficients $\{ \mu[k], \eta[k] \} = \{ a[k]+jb[k], c[k]+jd[k] \}$ are determined by gradient descent using

$$a^{(n+1)}[k] = a^{(n)}[k] + \varepsilon \frac {\partial C(\mu^{(n... ...(n)}[k] + \varepsilon \frac {\partial C(\mu^{(n)},\eta^{(n)})}{\partial b[k]} ,$$ (12)

$$c^{(n+1)}[k] = c^{(n)}[k] + \varepsilon \frac {\partial C(\mu^{(n... ...^{(n)}[k] + \varepsilon \frac {\partial C(\mu^{(n)},\eta^{(n)})}{\partial d[k]}$$ (13)

where $\varepsilon$ is the step size. The superscript $(n)$ corresponds to the value of $\mu$ and $\eta$ at the $n^{\text{th}}$ iteration. Each equation above represents a $3 \times 1$ vector equation³, i.e., one equation for each component of the vectors $a$, $b$, $c$, and $d$. The partial derivatives⁴ of $C(\mu,\eta)$ used in Eqs. 13 and 14 are

$$\frac {\partial C(\mu,\eta)}{\partial a[k]} = Re\{ I_0[k] + I_2[k]\}, \quad \frac {\partial C(\mu,\eta)}{\partial b[k]} = Im\{-I_0[k] + I_2[k]\},$$

$$\frac {\partial C(\mu,\eta)}{\partial c[k]} = Re\{ I_1[k] + I_3[k]\}, \quad \frac {\partial C(\mu,\eta)}{\partial d[k]} = Im\{-I_1[k] + I_3[k]\}$$

where

$$I_0[k] = \frac 4 {N_1 N_2 N_3} \sum_{n \in \Omega_d} \Big( \sigma \big( ... ...h_d[n,r]} + \chi\big(u_d[n,r] - \tilde{w}_d[n,r] \big) \Big) e^{j<x,\theta[k]>}$$

$$I_1[k] = \frac 4 {N_1 N_2 N_3} \sum_{n \in \Omega_d} \Big( \sigma \big( ... ...g_d[n,r]} + \chi\big(w_d[n,r] - \tilde{u}_d[n,r] \big) \Big) e^{j<x,\theta[k]>}$$

$$I_2[k] = 4\rho D^2[k] \mu[k]$$

$$I_3[k] = 4\rho D^2 [k] \eta[k].$$

The terms $I_0$ and $I_1$ are $3 \times 1$ vectors and are efficiently computed for all $k$ using 3D FFTs; each component of $I_0$ and $I_1$ requires a separate 3D FFT for a total of 6 FFTs. The notation $\nabla T_d \Big\vert _{N h_d[n,r]}$ represents the gradient of $T_d$ evaluated at the point $N h_d[n,r]$. The gradient is computed using symmetric difference equations⁵and trilinear interpolation is used to evaluate the gradient at the point $N h_d[n,r]$. The notation $\nabla S_d \Big\vert _{N g_d[n,r]}$ is defined similarly. The gradient descent insures that $\mu$ and $\eta$ have complex conjugate symmetry. Notice that $I_0$ and $I_1$ are the Fourier Transforms of a real signal implying that $I_0$ and $I_1$ have complex conjugate symmetry. This in turn implies that the gradients in Eqs. 13 and 14 have conjugate symmetry and therefore the basis coefficients $\mu$ and $\eta$ have complex conjugate symmetry.

The steps involved in estimating the basis coefficients $\mu$ and $\eta$ are summarized in the following algorithm.

Algorithm

1. Set $\mu[k] = \eta[k] = 0$ for $k \in \Omega_d$ and set $r = [1, 1, 1]^T$.
2. Compute $g_d^{-1}[n,r]$ using the procedure described in Section 2.3 and set $\tilde{w}_d[n,r] = g_d^{-1}[n,r] - \frac n N$.
3. Compute new forward basis coefficients $\mu[k] = a[k] + jb[k]$ using Eq. 13.
4. Compute the forward displacement field $u_d[n,r]$ using Eq. 11.
5. Compute $h_d^{-1}[n,r]$ using the procedure described in Section 2.3 and set $\tilde{u}_d [n,r] = h_d^{-1}[n,r] - \frac n N$.
6. Compute new reverse basis coefficients $\eta[k] = c[k] + j d[k]$ using Eq. 14.
7. Compute the reverse displacement field $w_d[n,r]$ using Eq. 11.
8. If the criteria is met to increase the number of basis functions then set $r = r +1$, and set the new coefficients in Eq. 11 to zero.
9. If the algorithm has not converged or reached the maximum number of iterations goto step 2.
10. Use forward displacement field $u_d[n]$ to transform the template image $T_d[n]$ and reverse displacement field $w_d[n]$ to transform the target image $S_d[n]$.

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