to resemble the wavelength of the laser (e.g., if the laser has color in the red spectrum, then consider only the red channel of the RGB image).
3.2 Pair-wise Optimization of Texture Registra-tion using Mutual InformaRegistra-tion
The estimation method of mutual information influences the accuracy and the complexity of registration, and may limit the possibilities of numerical optimiza-tion. If we refer to (3.2), the complexity comes from the term H(U, V) and any change in the transformation P modifies the joint r.v.
The entropy of a signal can be estimated using its histogram, which is a dis-crete approximation of the probability density function (pdf). Many techniques for estimation of mutual information rely on the joint histogram of the signals.
A simple estimation of the entropy can be obtained from the number of bins that are occupied.
The complexity of the estimation depends on how the joint histogram is pop-ulated. The first problem is to determine the entries in the joint histogram; if U and V are each sampled on a discrete grid, under the transformation T the grids may not be aligned; in the previous section, we considered signals defined on continuous domains, so we did not discuss this problem. The second problem is the quantization noise and other noise in the signals. From the perspective of [19] there are four general approaches to populate the joint histogram, two of them relevant for the first problem (grid alignment) and two for the second one (quantization and noise):
Standard Sampling (STD). The histogram is computed in the usual way, each sample adding a unit update for one histogram bin, using an interpo-lation method to estimate the signal on the second grid.
Partial Volume Estimation (PVE). Fractional updates are added to all bins corresponding to spatial neighbors of each sample (immediate neighbors, e.g., four, for a two-dimensional image), with the weights summing to 1 (see Partial Volume Interpolation (PVI) [47]). The method was generalized
3.2 Pair-wise Optimization of Texture Registration using Mutual Information
to PVE in [50], where the spatial neighborhood is selected using a general window.
Post-Parzen Windowing (PPZ). The histogram is constructed by standard sampling and then it is blurred (smoothed) with a convolution kernel.
In-Parzen Windowing (IPZ). Each sample (without being quantized) is con-voluted with the kernel and then the entries of the histogram are updated.
STD and PPZ do not achieve smooth estimations of the joint entropy (and mutual information). The PVE and IPZ achieve smooth estimations, with the degree of smoothness depending on the convolution kernel. If the kernel of PPZ and IPZ is the hat function with the support less or equal to the quantization size, they degenerate to STD. Similarly, if the spatial kernel of PVE is the hat function with support less than half of the spatial sampling, the sampling becomes nearest neighbor interpolation and PVE degenerates to STD.
PVE differences itself from the other methods in the sense that it does not rely on signal interpolation, although the grids are not aligned. The histogram obtained in this way could still be convoluted with a kernel to account for quan-tization and initial noise, therefore the four classes are not disjoint.
The complexity of the four methods is at least linear (STD) and at most quadratic in the number of samples N, depending on the support of kernels. If the support of the kernel is infinite (such as the Gaussian kernel) the complexity isN2; the sampling kernel (spatial) used in the PVI [47] determines also orderN complexity, with a constant factor depending on the dimensionality of the r.v.;
PPZ is cheaper than IPZ.
If the gradient and the Hessian of mutual information exist, they can be estimated with the same complexity, see [19] for the derivations.
The mutual information is known to be an objective function difficult to op-timize due to its multiple local maxima. Most nonlinear optimization techniques were applied to mutual information, as reported in literature (again, we refer to [56] for a synthetic review).
Several aspects should be considered when choosing an optimizer. It is impor-tant to note that the methods presented above are only estimations the mutual
3.2 Pair-wise Optimization of Texture Registration using Mutual Information
information. The quality of the estimate is better for more complex methods (such as IPZ). However, even if the data size is moderate, it is prohibitive to use a quadratic complexity algorithm (e.g., IPZ with a kernel with infinite support).
The conclusions of the evaluation presented in [19] indicate that IPZ with small and non-adaptive kernel performed worse than PVE and even STD.
Considering the increased complexity, multiresolution optimization is useful.
A special optimizer, similar to Levenberg-Marquardt but adapted to a multires-olution method using B-Spline kernels, was developed in [71]. The convergence rate is increased by the multiresolution strategy, and given good initialization and ideal conditions, it is claimed to have accuracy on the order of one hundredth of a pixel. Though, it is not easy to adapt multiresolution strategies to irregularly sampled signals, such as 3D surfaces.
3.2.1 Fast Gradient Estimation
In this thesis we used the optimization method proposed by Viola and Wells [79]. Although it is not based on the joint histogram, it can be considered an IPZ method and can be easily applied to irregularly sampled signals. It is based on stochastic sampling, and its time complexity, independent of the model’s or image’s resolution, is not prohibitive.
The core of the method is based on an approach to evaluate the entropy called EMMA [77], summarized by the following equations:
H(z) = −E(log(p(z))
In the first step, the entropy of an arbitrary r.v. zis estimated with a technique called Parzen windowing [20]. Given a finite samplingAofz, the pdf is estimated as a superposition of kernel functions G centered on the elements of A. In (3.4), the kernel G is parameterized by ψ; a common choice is the Gaussian kernel, in which case ψ is the associated covariance matrix. In the second step, the expectation of the logarithm is approximated with the average over a second sample set B.
3.2 Pair-wise Optimization of Texture Registration using Mutual Information
From (3.1), it turns out that in order to estimate the mutual information, we need to apply (3.4) for three r.v.: U, V and (U, V).
The advantage of this estimation method is that the entropy is theoretically differentiable with respect to z, the order depending only on the kernel function G. If the Gaussian kernel is chosen, then the entropy (and mutual information) is infinitely differentiable. A disadvantage of (3.4) is that its complexity is quadratic in the sampling size (NA·NB). To handle that, either the sampling sizes should be relatively small, or the kernel support should be bounded.
The fact that the r.v. V and (U, V) depend on P can be used to optimize P; indeed, changingP changes the U-V correspondence as well. Generically, the gradient of the entropy of a r.v. z(P) with respect to P can be found applying the differentiation chain rule to (3.2); if the kernel in the Parzen window is a Gaussian, then: Returning to the mutual information between U and V(P) from (3.1), its gradient with respect to P is:
d
dPI(U, V(P)) = dPd H(U) + dPd H(V(P))− dPd H(U, V(P))
= dPd H(V(P))− dPd H(U, V(P)) (3.6) In (3.6), one term was reduced, since it was independent of P. The most relevant term is the joint r.v. (U, V), which may change significantly even for small changes of P. The r.v. U is a tri-variate r.v. (the canonical form of the surface normal representation), and V is univariate (the intensity of an image).
Each sample of the r.v. (zi and zj from the sets A and B from (3.5)) con-tributes to the gradient. For the joint normal-intensity value (z = (U, V)), only the subvector corresponding toV, from (U, V), depends onP; (3.5) can be further expanded using:
V(P(X, Y, Z)) = V(x(P(X, Y, Z)), y(P(X, Y, Z)))
d
dPV = dxdV dPd x+ dydV dPd y (3.7)
3.3 Global Optimization of Texture Registration using Mutual