Having the grid and the deformation model, a cost function is designed which con-sists primarily of a similarity metric evaluating the difference between the source and the target images. The registration is then performed by minimizing the cost func-tion, which corresponds to minimizing the image dissimilarity between the deformed source image and the reference target image. For this purpose, different optimization algorithms can be used (Kleinet al., 2007) such as the Gauss-Newton method ( Bar-toli and Zisserman, 2004), the Levenberg-Marquardt algorithm (Kabus et al., 2004), or the method of conjugated gradients (CG) (Chun and Fessler, 2008;
Barber and Hose, 2005). During iterative optimization, the parameters of the model function are altered until the relative improvement becomes negligible.
Since in many cases the deformation is constrained, a regularization term is added to the cost function, weighted by a factor that counteracts the deformation.
Usually the regularization terms are designed such that the value of this term in-creases with the amplitude of the computed displacement. In the work of Chun and Fessler (2008, 2009b), a regularization method is presented based on a novel and very attractive penalization to ensure local invertibility. Alternatively, a rigidity penalty term can be used (Staringet al., 2007;ChunandFessler, 2009a). With regularization, the deformation is encouraged to be diffeomorphic and invertible (Ashburner, 2007).
A zoo of other registration approaches are existing. Template propagation for es-timating respiratory motion (R¨oschet al., 2002). The popular SIFT image features have been used to select control points (Franzet al., 2006). The Radon transform was also used for scale and rotation invariant image matching (Jiangsheng et al., 1998). The method of demons is also popular to track borders (Thirion, 1998). In this work, only two alternative implementations based on the optical flow principles have been considered.
Optical flow-based methods track the movement of image elements, assuming that each pixel is conserving its brightness value (Horn and Schunck, 1981; Ku-maret al., 1996;Broxet al., 2004). An alternative technique based on the popular block-matching scheme has also been implemented (Jain and Jain, 1981; Malsch et al., 2006;Chen, 2009). Block matching procedures are known to be fast and ex-periments with this specific implementation demonstrated very surprising accuracy.
The remainder of this chapter is structured as follows. A state of the art free-form deformation based registration technique is presented in section 8.2. This algorithm has been used for motion estimation in the experiments from the previous chapter.
For comparison a direct registration method based on block matching is described in section 8.3. Section 8.4 presents a comparative analysis of results obtained by the two techniques. An extensive discussion of parameters is provided. Further research and experiments are suggested in section 8.5.
The basis assumption of optical flow methods states that the brightnessE of one picture element does not change over time, hence
dE
dt = 0. (8.1)
Applying the chain rule for derivatives gives the following relation between infinites-imal spatial displacements and the change over time ofE:
∂Ex
∂x u+∂Ey
∂y v+∂Et
∂t = 0, (8.2)
withu= dxdt andv= dydt being the coordinates of the optical flow velocity vectors for the image element positioned in (x, y). Thus if the pixel intensityE does not remain constant over time, optical flow velocities can be computed according to the gradients in horizontal and vertical directions. Often image gradients are approximated by finite difference.
Based on the seminal paper of Horn and Schunk, various image registration methods were developed. One of these methods is the Sheffield image registration toolkit (ShIRT) from Barberand Hose (2005); Barber et al.(2007) on which a closer look was taken due to high expectations regarding its speed and robustness.
The following description of the registration method used in ShIRT borrows nota-tions from the original user manual. The implementation of the method is available on request to the author and has been developed since a decade.
The algorithm aims at minimizing a cost functionQ proportional to the dissim-ilarity between the images f and m = Γ (g), the source image g, deformed by the motion field Γ. The objective function can be written as
Q=X
(f−m)2. (8.3)
For each given pixel, it is assumed that any difference of intensities between f and m is due to an infinitesimal displacement of the pixel’s position. According to this fundamental principle of optical flow, a linear relashionship between the image difference and the horizontal and vertical displacements ∆x and ∆y can then be written as
f −m= ∆x 2
∂f
∂x+∂m
∂x
+∆y 2
∂f
∂y + ∂m
∂y
. (8.4)
The sum of squared differences used in Qis a simple and appropriate metric for capturing dissimilarity between binary images but tends to emphasize small relative differences between bright pixels. To overcome this problem, both source and target image are converted to binary imagesfb and mb, respectively, prior to registration.
The conversion to binary images is done by introducing a supplemental dimension for the intensity, effectively casting a 2D registration problem between a pair of planar images to a registration of volumetric 3D binary images:
fb−mb= ∆x 2
∂fb
∂x +∂mb
∂x
+∆y 2
∂fb
∂y +∂mb
∂y
+∆s 2
∂fb
∂s +∂mb
∂s
. (8.5)
The new dimension for the intensity channel is denoted bys. Since the intensity channel depends on the specific range of pixel’s values, it is appreciable to make this dimension independent from the spatial channels x and y. Therefore the difference between intensity channels is modelled as ∆s=s∆α+ ∆β.
The difference for continuous intensity images can now be expressed by integrat-ing oversthe difference between binary images, giving
f −m= ∆x 2
∂f
∂x+∂m
∂x
+∆y 2
∂f
∂y + ∂m
∂y
−∆α
2 (f +m)−∆β. (8.6) The objective is now to estimate smooth functions to represent the displacements
∆x, ∆y, ∆α, and ∆β. Thus the registration of planar 2D images boils down to the estimation of a smooth 4D function. To represent the deformation, a sparse grid of control points is introduced. For each control point located at the grid vertex (xi, yi), the bilinear kernel
φi(x, y) = max
0,1−
x−xi
d
×max
0,1−
y−yi
d
. (8.7)
is used to sample image elements at any given position (x, y).
The radius of the kernel depends on the spacing of control knots. The triangular basis function corresponds to the second order B-spline, many authors do prefer re-lying on cubic B-spline basis functions for their attractive smoothness and continuity properties. However, to the author’s knowledge, the practical relative advantage of cubic over linear B-splines has not yet been evaluated for image registration purpose.
Using such a grid for representing the deformation, the relation between image differences and the displacements of control points can be expressed by
f−m = X
i
axiφi
∂f
∂x +∂m
∂x
+ X
i
ayiφi
∂f
∂y + ∂m
∂y
− X
i
aαiφi(f−m)−2X
i
aαβiφi, (8.8) where the four scalar parametersaxi,ayi,aαi, andaβi must be estimated for every control pointi.
It is obvious that the number of control points influences the speed of the algo-rithm as well as the local accuracy of the displacements. If the number of control points is low, there are less computations to be executed. If on the other hand the spacing of these points is small, the displacements estimates becomes more accurate and even very slight motion can be captured.
Using vector and matrix notations, the relation can be written in shorten form:
f −m=aTT, (8.9)
whereaT represents the transposed vector of parameters to estimate and the matrix T represents the multiplication of basis functionsφi with sums of partial derivatives.
-15 -10 -5 0 5 10 15 -15-10-5 0 5 10 15 0
0.2 0.4 0.6 0.8 1
-15 -10 -5 0 5 10 15 -15-10-5 0 5 10 15 0
0.2 0.4 0.6 0.8 1
-15 -10 -5 0 5 10 15 -15-10-5 0 5 10 15 0
0.2 0.4 0.6 0.8 1
(a)falloff = 10−4 (b)falloff = 10−3 (c) falloff = 10−2 Figure 8.1: Various weightings of the image area surrounding the search block. A larger falloff factors of an exponential distribution allows selecting between uniform con-stant weighting where every pixel are equally important and point weighting where only the center pixel contributes to the evaluation of local image dissimilarity.
To enforce smoothness, hence, regularity of the deformation, a weighted regular-ization term is added which includes the square of the Laplacians. This regulariza-tion scheme was already proposed in the seminal paper of Horn and Schunk.
The regularized cost function can be expressed by Q=X
(f −m)2+λaTLTLa, (8.10) whereLstands for a discretized Laplacian operator andλis a user-defined parameter driving the strength of the regularization term.
In the implementation of ShIRT, the registration problem is iteratively solved by a conjugate gradient (CG) method. The current estimate for the parameters a at iterationt≥1 is noted byat. The accuracy of the estimate is iteratively updated asat+1 =at+ ∆at with
∆at=
TTT +λLTL−1
TT(f −mt)−λLTLat
. (8.11)