Distortion Mapping – An SVD HowTo

In this tech-talk, I would like to discuss how (and probably how not to) develop an optical distortion correction. It is a common optical problem, and one that actually can be compensated using an algorithm because, even if the image is distorted, it can still free from wavefront aberrations. I guess the title spills the beans already, it will be a follow-up on my love story with the singular value decomposition.

This problem is interesting, not only because it is common but also because it so clearly will demonstrate what the SVD is and its innate strengths. Now, before we dive in with some numbers, let’s just set some context.

Assumptions about the problem to solve

For a well aligned system with a large designed-in distortion, we have rotational symmetry, and radial polynomial will be the best way to address the effect. This is however not the problem I am interested in. The problem that I’d like to address is one where distortion is small and when we measure it, we cannot ignore the measurement noise. However, even if this distortion is small, it’s not small enough for us to ignore. Where would we find such an application? I think lithography may be one example where we may aim to correct it to better than one part in one hundred thousand. At this level of accuracy, we cannot assume anything about the nature of the distortion. The optics is as good as can be and there is still a distortion we can measure.

The distortion problem correction problem fits well into the SVD framework because we can collect data to fit the matrix formalism. Therefore, I will assume we will collect the distortion information on a grid that uniformly covers a rectangular image.

So let’s jump right in, both feet. Let’s look what we can expect from this approach. To generate some distortion data we can experiment with, I wrote a Matlab/Octave script. Readers that want to experiment with this can download it here.

N = 6;
[x,y] = meshgrid(-10:10, -13:13);
x = x/max(x(:));
y = y/max(y(:));
%a = randn(1,(N+1)*(N+2)/2).*(1:((N+1)*(N+2)/2))/N^2;
fx =  a(1)*x + a(2)*y + ...
      a(3)*x.^2 +  a(4)*x.*y + a(5)*y.^2 + ...
      a(6)*x.^3 +  a(7)*x.^2.*y + a(8)*x.*y.^2 + a(9)*y.^3 + ...
      a(10)*x.^4 + a(11)*x.^3.*y + a(12)*x.^2.*y.^2 + a(13)*x.*y.^3 + a(14)*y.^4 + ...
      a(15)*x.^5 + a(16)*x.^4.*y + a(17)*x.^3.*y.^2  + a(18)*x.^2.*y.^3 + a(19)*x.*y.^4 + a(20)*y.^5 ;
%     a(20)*x.^6 + a(21)*x.^5.*y + a(22)*x.^4.*y.^2 + a(23)*x.^3.*y.^3 + a(24)*x.^2.*y.^4 + a(25)*x.*y.^5 + a(26)*y.^6;

%b = randn(1,(N+1)*(N+2)/2).*(1:((N+1)*(N+2)/2))/N^2;
fy = b(1)*x + b(2)*y + ...
     b(3)*x.^2 + b(4)*x.*y + b(5)*y.^2 + ...
     b(6)*x.^3 + b(7)*x.^2.*y + b(8)*x.*y.^2 + b(9)*y.^3 + ...
     b(10)*x.^4 + b(11)*x.^3.*y + b(12)*x.^2.*y.^2 + b(13)*x.*y.^3 + b(14)*y.^4 + ...
     b(15)*x.^5 + b(16)*x.^4.*y + b(17)*x.^3.*y.^2 + b(18)*x.^2.*y.^3 + b(19)*x.*y.^4 + b(20)*y.^5;
%     b(21)*x.^6 + b(22)*x.^5.*y + b(23)*x.^4.*y.^2 + b(24)*x.^3.*y.^3 + b(25)*x.^2.*y.^4 + b(26)*x.*y.^5 + a(27)*y.^6;

fx = fx + 1*randn(size(fx))/400;
fy = fy + 1*randn(size(fy))/400;
[ux,sx,vx] = svd(fx);
[uy,sy,vy] = svd(fy);

SVD and Separation of Variables – An Analogy

The code above is just a MacLaurin expansion up to 5th order with random coefficients. The expansion defines displacements along the coordinate axis over a grid which we analyze using the singular value decomposition. The result of an SVD analysis is a matrix decomposition that splits the matrix into three terms,

The matrix “S” is diagonal and (conventionally) sorted in descending order. The columns of the U and V matrices are orthogonal and (usually) normalized,

Before we dive into the details of analyzing our model distortion, I would like to connect the SVD to something most readers may already be more familiar with, the method of separation of variables. When we want to solve a partial differential equation analytically, we assume that the solution can be written as,

This later leads us to a solution that we write as,

But this is exactly what the singular value decomposition offers us, a way to write our matrix M as a product of orthogonal functions independently spanning the row and column space of the matrix. When we prepared our inputs, we did something similar where our “basis functions” were all kinds of combinations of powers of x and powers of y multiplied by some (random) coefficient. The SVD analysis reduced our set of 20 input functions to two sets (U and V) of 6 functions. We can see this by running the script (a complete version of the script can be downloaded here)

The above results were generated without adding noise just to cleanly separate what is a solution and what is (numerically) zero.

Analyzing the distortion with noise

When we are faced with real data, at least for the assumptions I’ve made above, we cannot assume to know the size of our expansion. This is something we have yet to figure out. Therefore, lets add some noise and see how it looks,

With noise added, the structure of the singular values looks different. Where in the first case, the singular values dropped 12 orders of magnitude after the initial 6, there is now what looks like a cross over between two straight lines. This is actually quite common when we add white noise. For our distortion model, white noise is probably a good assumption because the measurements should be independent. The script now also shows us that our recovery of the distortion (up-sampled 32×32 times) is not perfect. The noise does not allow us to reliably recover more than 5 components. We could also argue, is the 5th component noise or is it data. By running the script enough times to gather some statistics, my own conclusion was that the error was on average lower if we consider it to be data and include it.

Comments on the script

There are some things in the script that may be worth pondering a little about. One that I find most disturbing is that when we interpolate and resample the left and right singular vectors (columns of the U & V matrices), they are no longer orthogonal. This means that we cannot use the interpolated vectors for projecting other distributions over this space. The lack of orthogonality does not appear to add any bias in the reconstruction though.

Also, the interpolations use the number of extracted singular values as the order of the polynomial used for interpolation. It is easy to argue for this choice for a model like this. This choice, however, would be one of the first I would try to confirm (or reject) when faced with real data.

If you found this interesting, please leave a comment. I wouldn’t mind rewriting this a an open-source C-program and include pre-warping of the input data. There are also some quite interesting numerical considerations how to structure the data to make an application of this theory able to handle mapping rates of giga pixels per second, which seems to be achievable with generic CPUs.

Thank you for reading.