The positions of the maxima of the electron density cannot be considered equal to the coordinates of the pixel with maximum electron density, because the resolu-tion of the grid is not sufficient (usually about 0.05 ˚A). An interpolation method must be used to locate the position of the maximum of the density more accu-rately. In addition to that, the t-sections of the superspace electron density have a general orientation with respect to the grid and their accurate extraction also demands interpolation of the density values between the grid points.
The basis of the interpolation method used inEDMA is a cubic spline inter-polation in one dimension (Press et al., 1996). The method has been generalized to arbitrary dimension using the philosophy of bicubic spline interpolation (Press et al., 1996). The principle of the method is illustrated in Fig. 2.4. The inter-polation of a n-dimensional electron density runs in n cycles. In every cycle, one-dimensional cubic spline interpolation is used to calculate a (k−1)D section of the kD density. This section contains the point, in which the density is to be determined. In the next cycle the newly calculated (k−1)D density is used as in-put and the dimension is further reduced. Ultimately, this leads the interpolated density of the point in question.
The result of the interpolation is independent of the order of axes, in which the density is interpolated. To prove this statement, we first construct a general formula for spline interpolation in one dimension. We have a function y(x) given by its tabulated values (xi, yi), i = 0. . . u −1, or, in short, by u-dimensional
38 CHAPTER 2. DEVELOPMENT OF BAYMEM
P[p , p ]1 2
Figure 2.4: Interpolation process in two dimensions. In first step, the rows of points xij, yij (circles) are interpolated along the dashed lines to obtain points of the vector
~yp1 (squares). This column is then interpolated along the full line to obtain the density at the point P (cross).
a set of n−2 linear equations of the form:
xk−xk−1
6 y00k−1+xk+1−xk−1
3 yk00+ xk+1−xk
6 y00k+1 = yk+1−yk
xk+1−xk − yk−yk−1
xk−xk−1
k = 1. . . u−2 (2.50) In addition to that, two more equations have to be supplied to obtain a unique solution. The most usual choice is y000 = y00u−1 = 0 to produce so called natural splines, but any other choice is possible. This set of linear equations can be compactly written in matrix form:
M~y00 =N~y
~y00=M−1N~y=S~y (2.51)
The matrixS is independent of the functional values yi and depends only on the values xi. In next step, we construct vectors A, ~~ B, ~C, ~D of dimension u from the functions A, B, C, D so that
A~·~y = Ayj
B~·~y = Byj+1
C~·~y = Cy00j D~·~y = Dyj+100
(2.52)
This is readily achieved by
Ak =
( 0 : k 6=j
A : k =j (2.53)
2.9. ANALYSIS OF ρM EM - PROGRAM EDMA 39 and analogously for B, ~~ C and D. With this notation, the interpolation formula~ can be written in a concise form as (compare Eqs. 2.48 and 2.51):
yp =A~·~y+B~·~y+C~TS~y+D~TS~y (2.54) or
yp =~s·~y (2.55)
The vector~s depends only on the values xi and xp.
Having obtained the general formula for the one-dimensional cubic splines we can turn to the higher-dimensional case. We will first explicitly consider the two dimensional case for the sake of simplicity. The generalization to the general nD case is then straightforward.
The 2D function is defined by its tabulated values xij,yij,i= 0. . . u−1, j = 0. . . v −1. We suppose that the values xij lie on a straight line, if one index is fixed and the other varies. This is equivalent to following condition: every four values xij, x(i+1)j, xi(j+1), x(i+1)(j+1) must form corners of a parallelogram.
(Fig. 2.4).
To interpolate a valueyp1p2 of the functiony(x1, x2) at point P[p1,p2], we first calculate a one-dimensional section of the two-dimensional density by interpola-tion along every row of pointsxij with the first index fixed. As a result, we obtain v pairs xpj, ypj, j = 0. . . v −1, where p denotes the first coordinate is equal to the coordinate p1 of the pointP. Using previous results (Eq. 2.55), we can write:
ypj =~s1·~yj, j = 0. . . v−1 (2.56) Notation ~yj denotes a column vector of dimensionu corresponding to the row of yij with j constant and i = 0. . . u−1. The symbol ~s1 stands for the vector ~s from Eq. 2.55 related to the rows of xij with the second index constant. Because of the condition on the regularity of the grid the vector~s1 is independent of the index j.
We form a one-dimensional vector~yp1 from the valuesypj andj Eqs. 2.56 can be written in a matrix form
~yp1 =Y~s1 (2.57)
where Y is the matrix of the values yij.
The final valueyp1p2 is obtained using Eq. 2.55 again:
yp1p2 =~s2T·~yp1 (2.58) Combining Eqs. 2.56 and 2.57, we obtain the final result:
yp1p2 =~s2TY~s1 (2.59) The proof follows immediately from Eq. 2.59. If the interpolation is inde-pendent of the order of the axes, the resulting formula must be invariant to the
40 CHAPTER 2. DEVELOPMENT OF BAYMEM transposition of the matrix Y and simultaneous exchange of the indices 1 and 2 of the vectors~s. Eq. 2.59 fulfills this requirement.
The matrix notation becomes impractical for more-dimensional cases. There-fore we turn to the tensor notation. Eq. 2.59 becomes (using the standard con-vention for summation over repeated indices):
yp1p2 =s1is2jyij (2.60) The general n-dimensional case follows immediately from the considerations in the previous paragraphs (and by analogy to Eq. 2.60):
yp1p2···pn =s1i1s2i2· · ·sninyi1i2···in (2.61) This equation is invariant to any permutation of the indices ik combined with corresponding permutation of the first indices of the numbers skik.
The interpolated value by the cubic-spline method is theoretically dependent on all values of the grid. In practice, the dependence strongly decreases with the distance from the central point. Therefore, in EDMA, every point is calculated from the pixels lying in some (user-definable) neighborhood of that pixel. This simplification does not introduce any detectable errors and saves computational time.