Alignment by Maximization of Mutual Information

“Alignment by Maximization of Mutual Information” by Paul Viola and William M. Wells III of the MIT Artificial Intelligence Laboratory was presented at the International Conference on Computer Vision (ICCV) in Cambridge, Massachusetts, in June 1995, and a journal version with the same title appeared in the International Journal of Computer Vision (IJCV) in 1997. It proposed aligning a model with an image by adjusting the transformation between them until the mutual information (MI) between model values and image intensities is as large as possible, so that nothing has to be known about how one is turned into the other: the lighting, the surface reflectance, or the physics of two different medical scanners. To make this practical, it estimated the entropies from small random samples with Parzen windows and climbed the MI by stochastic gradient ascent. The paper demonstrated the method on magnetic resonance (MR) images, on 3D object models aligned to cluttered video, on head tracking and on a view-based model. With contemporaneous work from Leuven, it is among the first uses of MI as a registration criterion; Zitová and Flusser’s (2003) survey calls it one of the first. The image registration article covers MI as a similarity measure as it is used today; this article covers the publication.

This article is read in the IJCV 1997 version, in the authors’ copy from Wells’s MIT page (a 29-page manuscript headed with the journal reference), and attributes method details, equation numbers and results to that version. The conference paper was read in the authors’ preprint marked “To appear in International Conference on Computer Vision 95”, archived by archive.org from Viola’s MIT AI Laboratory page; it differs from the journal version in ways listed under Versions, and it may also differ from the printed proceedings, which we could not read.

Problem

The journal paper frames alignment with an imaging equation. A model $u(x)$ and an image $v(y)$ are related by

$$v(T(x)) = F(u(x), q) + \eta,$$

where $T$ is the transformation, or pose, that maps model coordinates into the image, $F$ is the imaging function, $q$ collects exogenous parameters such as the lighting, and $\eta$ is noise, usually assumed Gaussian. Alignment is hard, the authors argue, because $F$ is difficult to model, $q$ is usually unknown, and the space of $T$ is large: six dimensions for a rigid object and in principle unbounded for a non-rigid one.

Their running example is a 3D model of a person’s head, with a surface normal at each point, and a video image of the head. Rendered under a Lambertian assumption and lit from the right, the model at the correct pose looks to a person very like the video image; yet the paper reports that the pixels of the real image and the correctly posed rendering are uncorrelated. The paper’s discussion contrasts two common ways of comparing a model with an image and rejects both:

• Correlation compares intensities directly. The paper notes that correlation matching assumes a Gaussian conditional distribution of the image given the model, and that correlation assumes model and image are identical or related by a linear function; it argues that negating the image or applying a non-linear intensity map defeats conventional correlation.
• Edge-based matching (edge detection followed by a comparison of edge maps), such as Huttenlocher’s Hausdorff distance, compares edges or gradient magnitudes. It works for objects with sharp creases or useful silhouettes, but assumes that the edges are stable under lighting changes and that the model is well described by edges. The paper’s face model, it notes, has no creases or points of high curvature from which such a representation could be built.

The same problem appears in medical imaging, which the paper names as a target: data from one sensor, such as MR, must be aligned with data from another, such as computed tomography (CT), whose intensities are not a function of each other. The journal paper’s acknowledgments say the authors were partly inspired by Hill and Hawkes’s work on medical registration, and its related-work section describes Hill, Studholme and Hawkes (1994) as aligning MR, CT and other modalities with third-order moments of the joint histogram.

Contribution

The journal version presents four things:

1. Mutual information as an alignment criterion, derived from maximum likelihood by successive relaxations: from a known imaging function to an unknown but consistent one, from consistency to conditional entropy, and from conditional entropy to MI. The criterion assumes only, in the paper’s words, that “the model tells more about the scene” when correctly aligned.
2. A differentiable entropy estimate from samples, which the authors call EMMA, for “EMpirical entropy Manipulation and Analysis”. It needs no parametric form for the distributions.
3. Stochastic maximization: each step uses a gradient estimated from a few dozen random samples, which keeps the cost per step small and, the authors argue, helps the search past small local maxima.
4. Experiments in four settings: MR–MR alignment, a 3D skull model aligned to cluttered video, head tracking in video, and alignment of a novel image to a model made of two images.

The paper describes the result as intensity-based rather than feature-based, and as more robust than correlation, because MI is unchanged if the image is negated and tolerates a wide range of non-linear intensity maps.

Method

From likelihood to mutual information

The derivation starts from the likelihood of the image given the model and pose, which in full requires integrating over all imaging functions and exogenous parameters. The paper sets this aside as probably infeasible. A nested search that alternately fits $F$ and $q$ and re-estimates $T$ is possible, but needs assumptions about the form of $F$ and a strong prior on its smoothness to converge. It replaces the search for $F$ with a test of consistency: under a correct pose, points with similar model values should have similar image values. Its first measure, a weighted sum of squared intensity differences between pairs of points that are close in model value (its eq. (7)), has an obvious flaw that the paper points out: it is maximized by constancy, for instance by a scale that maps the whole model onto a single image point. A weighted neighbor approximation of $F$, a kernel regression of image values on model values (eq. (8)), gives a likelihood that clearly separates the aligned head data from data with the model 3 cm too low: log likelihoods of 1079.49 against 537.34.

The paper then lists three drawbacks of weighted neighbor likelihood: it works only when the image is a function of the model, which, as the paper notes, CT and MR are not of each other; under a Gaussian noise model, occlusion and specular highlights ruin otherwise good matches; and it still favors constant solutions. The paper argues that moving to entropy removes the first two. Under the condition that the conditional distribution of the image matches the assumed noise, the conditional log likelihood equals $-N$ times the conditional entropy of the image given the model, for $N$ model points (eq. (11)); estimating the conditional entropy directly drops that condition. Conditional entropy is still low when the image region is itself constant, so the paper maximizes the mutual information

$$I\big(u(x), v(T(x))\big) = h\big(u(x)\big) + h\big(v(T(x))\big) - h\big(u(x), v(T(x))\big),$$

its eq. (13), with $h$ the differential entropy and $x$ a random point of the model. The first term does not depend on $T$. The second rewards poses that project the model onto complex, high-entropy parts of the image, and the third, the negative joint entropy, rewards poses under which the model explains the image; together, in the authors’ summary, they “find complexity and explain it well”. The second term is what stops the collapse onto a uniform region that defeated consistency.

Entropy from samples

The estimate replaces the unknown density with a Parzen window estimate, a sum of Gaussian kernels $G_\psi$ with covariance $\psi$ centred on the points of a random sample $A$, and replaces the expectation in $h(z) = -E[\ln p(z)]$ with an average over a second sample $B$ (eq. (15)):

$$h(z) \approx -\frac{1}{N_B} \sum_{z_i \in B} \ln \frac{1}{N_A} \sum_{z_j \in A} G_\psi(z_i - z_j).$$

Any differentiable kernel would do; the paper suggests the Cauchy density as another choice. Differentiating with respect to the pose gives (eq. (16))

$$\frac{d}{dT} h(z) \approx \frac{1}{N_B} \sum_{z_i \in B} \sum_{z_j \in A} W_z(z_i, z_j)\, (z_i - z_j)^{\top} \psi^{-1} \frac{d}{dT}(z_i - z_j), \qquad W_z(z_i, z_j) = \frac{G_\psi(z_i - z_j)}{\sum_{z_k \in A} G_\psi(z_i - z_k)}.$$

The weight $W_z$ is a softmax over negative Mahalanobis distances, which the paper relates to the softmax of neural networks: it is near 1 when $z_i$ is much closer to $z_j$ than to any other element of $A$. Reducing entropy therefore means pulling together values that are already near each other in value, tightening the clusters. With block-diagonal kernel covariances for the joint density, the MI gradient (eq. (17)) has two parts: one pushes apart samples that are close in image intensity, the other pulls together samples that are close in both image intensity and model value. Distances are measured between values, not image positions; for a linear $T$, $\frac{d}{dT} v(T(x_i)) = \nabla v(T(x_i))\, x_i^{\top}$, so the pose gradient comes from the image gradient at the projected sample points.

Stochastic maximization

The algorithm repeats three steps: draw a sample $A$ of size $N_A$, draw a sample $B$ of size $N_B$, and update $T \leftarrow T + \lambda\, \widehat{dI/dT}$ with learning rate $\lambda$. Using every model point would cost time quadratic in the number of points; small samples give noisy but cheap gradient estimates, an instance of stochastic approximation, for which the paper cites Widrow and Hoff (1960) and textbooks by Ljung and Söderström and by Haykin. The paper gives $N_A = N_B = 50$ as an example of a sample size that suffices in practice, states that the authors have proved convergence to a pose close to a local optimum (referring to Viola’s thesis), and reports that the sampling noise helps the search escape small local minima. The kernel covariances, assumed diagonal, are chosen to maximize the likelihood of the samples, which is equivalent to minimizing the cross entropy with the true distribution, and are estimated on-line by a similar stochastic scheme. Bilinear interpolation is used wherever the image is indexed at non-integer positions.

Worked Example: Stochastic Alignment Under a Non-Monotonic Intensity Map

The example below implements eqs. (15)–(17) and the stochastic update in NumPy and SciPy on a synthetic $128 \times 128$ image. The second image is the first moved by a rigid pose (a 6° rotation and a translation) and then passed through a quadratic intensity map, whose centre is chosen so that the two images are uncorrelated at the true pose. Each run starts from a random pose and takes 500 steps at a learning rate of 0.003 and then 500 at 0.0003. The two-stage schedule follows the paper’s MRI experiment; the rates are ours, because they depend on how the pose is parameterized. The sample size $N$ is varied.

import numpy as np
from scipy import ndimage

# Model image u: smooth blobs and a square on a 128 x 128 grid, intensities in [0, 1]
n = 128
yy, xx = np.mgrid[0:n, 0:n] / (n - 1) * 2 - 1            # pixel centres in [-1, 1]
u_img = np.zeros((n, n))
for cx, cy, s, a in [(-.4, -.3, .25, 1.0), (.35, .1, .18, .7), (0, .5, .3, .5)]:
    u_img += a * np.exp(-((xx - cx)**2 + (yy - cy)**2) / (2 * s**2))
u_img[(abs(xx - .3) < .25) & (abs(yy + .45) < .2)] += .6
u_img = ndimage.gaussian_filter(u_img, 2)
u_img = (u_img - u_img.min()) / np.ptp(u_img)

def sample(img, pts):                    # bilinear lookup at [-1, 1] coordinates
    c = (pts + 1) / 2 * (n - 1)
    return ndimage.map_coordinates(img, [c[:, 1], c[:, 0]], order=1, mode="nearest")

def pose(p):                             # p = (angle, tx, ty) -> 2 x 3 matrix T
    c, s = np.cos(p[0]), np.sin(p[0])
    return np.array([[c, -s, p[1]], [s, c, p[2]]])

# Image v: u moved by a rigid pose, then a non-monotonic (quadratic) intensity map
# whose centre m is chosen so that v is uncorrelated with u at the true pose
p_true = np.array([np.deg2rad(6), 0.08, -0.06])
T = np.vstack([pose(p_true), [0, 0, 1]])
grid = np.stack([xx.ravel(), yy.ravel(), np.ones(n * n)])
u_moved = sample(u_img, np.linalg.solve(T, grid)[:2].T).reshape(n, n)
ms = np.linspace(0, 1, 1001)
m = ms[np.argmin([abs(np.corrcoef(u_img.ravel(), (u_img.ravel() - mm)**2)[0, 1]) for mm in ms])]
v_img = (u_moved - m)**2 / max(m, 1 - m)**2
gy, gx = np.gradient(v_img, 2 / (n - 1))                  # image gradient in [-1, 1] units

def mi_gradient(p, rng, N, sigma=0.1):
    """Stochastic estimate of dI/dp, eq. (17) of the IJCV 1997 version."""
    xa, xb = rng.uniform(-.8, .8, (N, 2)), rng.uniform(-.8, .8, (N, 2))   # samples A, B
    c, s = np.cos(p[0]), np.sin(p[0])
    def values(x):
        y = np.c_[x, np.ones(N)] @ pose(p).T
        dy = np.stack([np.c_[-s * x[:, 0] - c * x[:, 1], c * x[:, 0] - s * x[:, 1]],
                       np.tile([1., 0], (N, 1)), np.tile([0., 1], (N, 1))], 1)
        return sample(u_img, x), sample(v_img, y), np.einsum('ik,ipk->ip',
                    np.c_[sample(gx, y), sample(gy, y)], dy)     # dv/dp per sample
    ua, va, ja = values(xa)
    ub, vb, jb = values(xb)
    dv, du = vb[:, None] - va[None, :], ub[:, None] - ua[None, :]   # i in B, j in A
    Wv = np.exp(-dv**2 / (2 * sigma**2)); Wv /= Wv.sum(1, keepdims=True)
    Wuv = np.exp(-(dv**2 + du**2) / (2 * sigma**2)); Wuv /= Wuv.sum(1, keepdims=True)
    k = dv * (Wv - Wuv) / sigma**2
    return (k.sum(1) @ jb - k.sum(0) @ ja) / N

def align(p, N, seed):
    rng = np.random.default_rng(seed)
    for rate, iters in [(0.003, 500), (0.0003, 500)]:     # two stages, as in the paper
        for _ in range(iters):
            p = p + rate * mi_gradient(p, rng, N)
    return p

def error_px(p):                        # mean displacement of the image corners, pixels
    c = np.array([[-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]])
    d = c @ (pose(p) - pose(p_true)).T
    return np.linalg.norm(d, axis=1).mean() * (n - 1) / 2

r = np.corrcoef(u_img.ravel(), sample(v_img, (T @ grid)[:2].T))[0, 1]
print(f"correlation of u and v at the true pose: {r:+.3f}")
for N in (10, 20, 50):
    start, final = [], []
    for trial in range(20):
        rng = np.random.default_rng(100 + trial)
        p0 = p_true + np.r_[np.deg2rad(rng.uniform(-15, 15)), rng.uniform(-.15, .15, 2)]
        start.append(error_px(p0))
        final.append(error_px(align(p0, N, seed=trial)))
    f = np.array(final)
    print(f"N = {N:2d}: start {min(start):4.1f}-{max(start):4.1f} px, final median "
          f"{np.median(f):5.2f} px, max {f.max():6.2f} px, within 1 px: {(f < 1).sum()}/20")

Output (NumPy 2.3.5, SciPy 1.16.3; two runs gave identical output):

correlation of u and v at the true pose: +0.004
N = 10: start  7.0-24.0 px, final median  0.95 px, max  31.08 px, within 1 px: 12/20
N = 20: start  7.0-24.0 px, final median  0.41 px, max   0.96 px, within 1 px: 20/20
N = 50: start  7.0-24.0 px, final median  0.32 px, max   0.50 px, within 1 px: 20/20

At the true pose the correlation between the two images is essentially zero, so a correlation score has nothing to climb, yet the MI gradient brings every one of the 20 starting poses to within a pixel with 20 or 50 samples per set. Ten samples per set make the gradient so noisy that some runs stop short of a pixel and one wanders far from the true pose; 50 samples tighten the final poses further. The cost of a step grows with $N_A N_B$, so the sample size trades the cost of each step against the noise in it, the trade-off the paper describes. The example uses a 3-parameter rigid pose to stay short, where the paper’s MRI experiment used a 6-parameter affine transform, and a single kernel width for the marginal and joint densities, as that experiment did.

Results

All timings in the journal version are for a Sun Sparc 10 workstation.

MR images

The first experiment aligned the proton-density and T2-weighted images of a double-echo MR scan, acquired simultaneously so that the correct alignment is close to the identity. The transform was a 6-parameter affine map; intensities were scaled to $[0, 1]$, the kernels had $\sigma = 0.1$, samples had 20 points, and the learning rate was 0.02 for 500 iterations and then 0.005 for 500, in 12 s in total. From 50 random initial poses that varied by 32 pixels in position (a little under a third of the head’s width), by 28° in rotation and by up to 20% in scale, the paper reports that the correct alignment was found reliably, with final alignments well within one pixel in position and within 0.5% of the identity matrix in rotation and scale.

The authors concede that normalized correlation would probably have aligned this pair too, since the two images are similar. They then passed the model image through a non-monotonic (quadratic) intensity transformation, which the conference preprint specifies as subtracting 0.5, squaring and rescaling to $[0, 1]$, and report that alignment performance was not significantly affected. Correlation, they note, would fail on this pair, since the two images are more similar in value when badly misaligned and not overlapping than when aligned.

3D skull model to video

A 3D skull model with surface normals, derived from a CT scan, was aligned to cluttered video images of the skull, without segmenting the images. The model values are two components of the surface normal and the image values are video brightness. The pose was a double quaternion followed by perspective projection for a camera with an 18° viewing angle. Kernel variances were 0.3 for each normal component, 0.2 for intensity and 0.15 for the image entropy estimate; samples had 50 points; hidden points were pruned with a Z-buffer every 200 iterations; and because a rotation step moves points about 100 times further than a translation step of the same size for a model with a radius of about 100 mm, rotation learning rates were set 100 times lower than translation rates (0.0005 and 0.05 for 2,000 iterations, then 0.0001 and 0.01 for 2,000 more), for about 30 s per alignment.

Table 1 of the journal version reports four experiments of 50 random initial poses each. From a base pose found by an initial alignment, each translation axis was offset by a uniformly distributed amount up to $\Delta T$ and the model was rotated about a random axis by up to $\Delta\theta$. The spread of the model centroid is given in the table as $\sigma_X$, $\sigma_Y$, $\sigma_Z$ in millimetres (the text calls these variances, but the symbols and units are those of standard deviations), the mean angular error $|\Delta\theta|$ in degrees, and the success rate; final statistics are computed over the successful runs only:

$\Delta T$ (±mm)	$\Delta\theta$ (°)	Initial $\sigma_X$, $\sigma_Y$, $\sigma_Z$	Initial $\lvert\Delta\theta\rvert$	Final $\sigma_X$, $\sigma_Y$, $\sigma_Z$	Final $\lvert\Delta\theta\rvert$	Success
10	10	5.94, 5.56, 6.11	5.11	0.61, 0.53, 5.49	3.22	100%
30	10	16.53, 18.00, 16.82	5.88	1.80, 0.81, 14.56	2.77	96%
20	20	10.12, 12.04, 10.77	11.56	1.11, 0.41, 9.18	3.31	96%
10 to 20	20 to 40	14.83, 15.46, 14.466	28.70	1.87, 2.22, 14.19	3.05	78%

The table shows something the text does not discuss: the spread in $X$ and $Y$ shrinks to one or two millimetres, but the spread in $Z$, the depth axis, hardly changes, and the mean angular error stays near 3°. By our reading, the video constrains the skull’s position across the image much better than its distance from the camera. The text says alignment is most sensitive to rotation in depth, because only visible points contribute to the gradient; that translations as large as half the skull’s diameter and in-plane rotations up to 45° can be accommodated; and that outside this capture range gradient ascent alone does not converge. The largest translation offsets in the table are 30 mm per axis, at most about 52 mm in all three axes together; half the diameter is about 100 mm if the skull has the radius of about 100 mm that the paper gives for its models (our arithmetic), so that claim rests on experiments the table does not report. With an artificial occlusion covering the chin, the recovered pose was very close to the correct one, though the authors note that occluded images can take longer to converge.

Head tracking

A 3D head model with normals from a Cyberware scan was aligned to each frame of a video sequence digitized at 3 frames per second, with the previous frame’s pose as the starting point and a manual starting pose for the first frame; each refinement took about 10 s. The authors list how this differs from the skull experiments: the face has no creases for an edge-based method to use, its albedo varies (eyebrows, lips), and the model is cut off at the chin and forehead, so the occluding contours of model and image do not agree. The results are shown as figures only.

View-based alignment

When no 3D model exists, the model can be a set of images. The paper shows two images of a car in the same pose under different lighting that are anti-correlated, which no variant of correlation can match, and a third, novel image whose scatter plot against either model image shows no consistent relationship at all, so that neither correlation nor EMMA can match it to one image. Its solution draws on photometric stereo: since two images under known lighting determine the surface normals, a pair of model images should carry information about any third image, even though the photometric stereo is never computed. Using the pair as a two-component model value, with $\sigma = 0.1$ and 1,000 iterations at a rate of 0.002, randomized experiments gave a capture range of about 40% of the car’s length and width and 35° of rotation.

Versions

• ICCV 1995, pages 16–23 of the proceedings, the paper that received the awards described under Impact. Its DOI leads to IEEE Xplore, and the open-access link that OpenAlex gives for it is Viola’s thesis. We read instead the authors’ 23-page preprint marked “To appear in International Conference on Computer Vision 95”, linked from Viola’s MIT page and captured by archive.org in November 1997. It starts directly from MI, without the journal’s derivation from likelihood, and writes $H$ for entropy; its entropy estimate, gradient and stochastic update are those of the journal’s eqs. (15)–(17). Its experiments are the same MR alignment with the same settings; a synthetic surface with three Gaussian bumps aligned to images rendered under different lighting (its Table 1: three experiments of 50 poses, with 100%, 100% and 80% success); the skull, whose Table 2 has the first three rows of the journal’s Table 1, in another order and without the fourth; and the car, with samples of 20 points. It has no head tracking. Where the journal reports a convergence proof, the preprint says the authors had been unable to prove that the derivative estimates are unbiased. Its experiment checklist gives the MR learning rates as 0.02 and then 0.0005, against 0.005 in its own MR section.
• Viola’s PhD thesis (MIT, submitted in June 1995 and issued as AI Technical Report 1548), with the same title, is the longer account to which the journal paper refers for its convergence proof and for more on robustness. We read its abstracts, in the copy on Viola’s archived Microsoft Research page. Unlike the journal version, its report abstract lists the registration of MR with CT images among the experiments, and its thesis abstract claims that stochastic approximation speeds up alignment by at least a factor of 500 over gradient descent; we did not check either in the body of the thesis.
• IJCV 1997, volume 24, issue 2, pages 137–154, read in the authors’ manuscript. Springer’s copy is not open access and could not be compared with it. Its abstract lists the MR experiments, 3D model alignment with clutter and occlusion, head tracking and view-based alignment.

The journal authors’ copy contains a sign error. Its eq. (12) writes the conditional entropy as $h(v \mid u) = h(u) - h(u, v)$; the correct form is $h(u, v) - h(u)$, so the printed expression is the negative of the conditional entropy. The derivation from likelihood, and with it this equation, is not in the conference preprint. The error does not reach the method, which maximizes MI in eq. (13). Eq. (10) also introduces the joint entropy of “x and y” while writing $h(z, y)$.

Impact

The PAMI Technical Committee lists the paper under the ICCV Best Paper Honorable Mention Award for 1995, one of three 1995 papers so listed. It also lists the paper among the eleven winners of the Helmholtz Prize for 2013, eighteen years after the conference; the committee describes the biennial prize, presented at each ICCV, as recognizing ICCV papers from ten years earlier with significant impact on computer vision research.

The paper’s main influence was on medical image registration. The authors applied the method to 3D volumes in Medical Image Analysis in 1996 (Wells et al.), registering MR with CT and with positron emission tomography (PET) images, and describing surgical uses. In Leuven, Maes, Collignon and colleagues (1997) state that their paper expands on ideas first presented by Collignon et al. at the Information Processing in Medical Imaging conference, which their reference list dates to June 1995, the same month as ICCV, and name Viola and Wells’s ICCV and Medical Image Analysis papers as related work. The journal paper, for its part, describes an earlier 1995 paper by Collignon et al., at the CVRMed conference, as using joint entropy as a criterion for CT–MR registration and probing it without describing a search method. Zitová and Flusser’s 2003 survey of image registration calls MI methods the leading technique in multimodal registration and names the IJCV paper as one of the first articles proposing the technique.

Limitations

• A local method. The authors state that alignment is reliable only from a starting pose near the correct one. In the skull experiments, success falls to 78% for the largest initial errors, and the head tracker relies on the previous frame’s pose.
• Depth and occlusion. Alignment is most sensitive to rotation in depth, because hidden points give no gradient, and the depth spread in Table 1 barely shrinks. Occluded images take longer to converge.
• Hand-set parameters. Kernel widths and learning rates differ between experiments, and the 3D experiments need separate rates for rotation and translation; the covariances can be estimated on-line, but the learning rates are chosen by hand.
• The relationship must exist in the data. MI needs a consistent statistical relationship between model and image values; the view-based experiment shows a novel image for which a single model image has none.
• Evaluation. The success rates in Table 1 are measured against a base pose found by the method itself, and the MR experiment against the identity transform expected from simultaneous acquisition. By our reading, this measures repeatability and capture range rather than accuracy against an independent reference; Maes et al. (1997), by contrast, validated their histogram-based MI against a stereotactic reference solution.

What Came After

The Leuven paper estimated MI from the joint histogram of the overlapping volumes rather than from Parzen windows, maximized it with Powell’s direction-set method, introduced partial-volume interpolation so that the histogram changes smoothly with the pose, and validated rigid registration of CT, MR and PET images against a stereotactic reference solution, reporting subvoxel accuracy. Its authors noted that comparing their histogram estimates with Viola and Wells’s Parzen windowing needed more experiments. Zitová and Flusser (2003) summarize Studholme, Hill and Hawkes (1999), whose title calls their measure overlap invariant, as comparing joint entropy, MI and normalized MI, the last a new idea claimed to outperform the others, with histogram estimates instead of the Parzen windows of earlier work; they describe Thévenaz and Unser as keeping Parzen window estimates but maximizing MI with the Jeeves or Levenberg–Marquardt methods and speeding it up with spline pyramids. The image registration article covers MI and normalized MI as registration criteria today, with their histogram estimates, interpolation and search.