Bilinear wavefront transformation

Truncated expansions such as Zernike polynomials provide a powerful approach for describing wavefront data. However, many simple calculations with data in this form can require significant computational effort. Important examples include recentering, renormalizing, and translating the wavefront data. This paper describes a technique whereby these operations and many others can be performed with a simple matrix approach using monomials. The technique may be applied to other expansions by reordering the data and applying transformations. The key is the use of the vectorization operator to convert data between vector and matrix descriptions. With this conversion, one-dimensional polynomial techniques may be employed to perform separable operations. Examples are also given for differentiation and integration of wavefronts.

K. Dillon, “Bilinear wavefront transformation,” J. Opt. Soc. Am. A, vol. 26, no. 8, pp. 1839–1846, 2009. (pdf)

For code, click more…

The Matlab code for the functions and matrices from the paper is given here, where “N” is the polynomial order.

Code to generate permutation matrices:

P = zeros(1/2*(N+1)*(N+2),(N+1)^2);
for i=0:N,
  for j=0:N-i,
    k = 1/2*(i+j)*(i+j+1)+j;
    k0 = i+(N+1)*j;
    P(k+1,k0+1)=1; 
  end;    
end;

Code to generate matrix for scaling by “alpha”:

Ds = diag(alpha.^(0:N));

Code to generate matrix for translation by “a”:

Dt = zeros(N+1, N+1);
for r=1:N+1,
 for c = 1:r,
  Dt(r,c) = a^(r-c)*factorial(r-1)/factorial(c-1)/factorial(r-c);
 end;
end;

Code to generate differentiation matrix:

Dd = diag((1:N), -1);

Code to generate integration matrix:

Di = diag((1:N).^(-1), 1);

To pad a coefficient vector (Zernike or monomial) up to the next higher order for integration, we simply append zeros to it as in the following:

c = [c; zeros(N+2,1)];

To form the vectorized transformation matrix, we first form the matrices for each coordinate. For example if we want to scale then translate by different amounts in both coordinates we would form:

Dx = Dt1*Ds1;
Dy = Dt2*Ds2;

With only a single-coordinate transformation, such as if we want to integrate in the x direction, we form:

Dx = Di;
Dy = eye(N+1);

and vice versa for the integration in y. Differentiation is treated similarly.
Finally the vectorized transformation matrix can be computed as:

Mm = P*kron(Dy.', Dx.')*P.';

Given a matrix T which converts monomials to Zernike coefficients, we may form the Zernike transformation matrix as:

Mz = inv(T)*P*kron(Dy.', Dx.')*P.'*T;

An algorithm for computing T is:

T = zeros(1/2*(N+1)*(N+2),1/2*(N+1)*(N+2));
for n=0:N,
  for m=-n:2:n,
    k_z=1/2*(n*(n+2)+m);
    d = m<0;
    if m==0, a=sqrt((2*n+2)/2); else a=sqrt(2*n+2); end;
      for l=0:(abs(m)-d)/2,
        for s=0:(n-abs(m))/2,
          for t=0:(n-abs(m))/2-s,
          b = (-1)^(s+l)*factorial(n-s)/factorial(s) ...
            /factorial((n+abs(m))/2-s)/factorial((n-abs(m))/2-s) ...
            *choose(abs(m), 2*l+d)*choose((n-abs(m))/2-s,t);
          i = n-d-2*(s+t+l);
          j = 2*l+d+2*t;
          k_m=1/2*(i+j)*(i+j+1)+j;
          T(k_m+1,k_z+1) = T(k_m+1,k_z+1) + a*b;
        end;
      end;
    end;
  end;
end;