Cartographie des distorsions - Un mode d'emploi de la SVD

Dans cet exposé technique, j'aimerais expliquer comment (et probablement comment ne pas) mettre au point une correction de la distorsion optique. Il s'agit d'une problème optique courantIl s'agit d'une aberration de front d'onde, qui peut être compensée par un algorithme, car même si l'image est déformée, elle peut toujours être exempte d'aberrations de front d'onde. Je suppose que le titre annonce déjà la couleur, il s'agira d'une suite à mon histoire d'amour avec la décomposition en valeurs singulières.

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.

Hypothèses sur le problème à résoudre

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.

Le problème de la correction des distorsions s'inscrit bien dans le cadre de la SVD, car nous pouvons collecter des données qui s'adaptent au formalisme matriciel. Par conséquent, je supposerai que nous collecterons les informations de distorsion sur une grille qui couvre uniformément une image rectangulaire.

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 who want to experiment with this can download it ici. (Vous pouvez envoyer un courriel à contact to get access if you don’t already have it)

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

Le code ci-dessus est simplement une expansion de MacLaurin jusqu'au 5e ordre avec des coefficients aléatoires. L'expansion définit les déplacements le long de l'axe des coordonnées sur une grille que nous analysons à l'aide de la décomposition en valeurs singulières. Le résultat d'une analyse SVD est une décomposition matricielle qui divise la matrice en trois termes,

The matrix “S” is diagonal and (conventionally) sorted in descending order. The columns of the U et V sont orthogonales et (généralement) normalisées,

Avant de nous plonger dans les détails de l'analyse de la distorsion de notre modèle, j'aimerais faire le lien entre l'UDS et quelque chose que la plupart des lecteurs connaissent peut-être déjà, la méthode de séparation des variables. Lorsque nous voulons résoudre une équation différentielle partielle de manière analytique, nous supposons que la solution peut être écrite sous la forme suivante,

Cela nous conduit ensuite à une solution que nous écrivons comme suit,

C'est exactement l'objectif de la décomposition en valeurs singulières. Elle permet d'écrire notre matrice M as a product of functions which independently span 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 et des pouvoirs de y multipliée par un coefficient (aléatoire). L'analyse SVD a réduit notre ensemble de 20 fonctions d'entrée à deux ensembles (U et V) de 6 fonctions. Nous pouvons le constater en exécutant le script (une version complète du script peut être téléchargée). ici)

Les résultats ci-dessus ont été générés sans ajouter de bruit afin de séparer clairement ce qui est une solution et ce qui est (numériquement) zéro.

Analyse de la distorsion avec le bruit

When we face 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,

Avec l'ajout de bruit, la structure des valeurs singulières semble différente. Alors que dans le premier cas, les valeurs singulières chutaient de 12 ordres de grandeur après les 6 initiales, on observe maintenant ce qui ressemble à un croisement entre deux lignes droites. Ce phénomène est en fait assez courant lorsque l'on ajoute du bruit blanc. Pour notre modèle de distorsion, le bruit blanc est probablement une bonne hypothèse car les mesures devraient être indépendantes.

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 lower if we consider it to be data and include it. Even more valuable is, the fact that by including that 5th point, the mapping suppressed outliers which would otherwise appear 1.5% or 2% of the times we rerun the simulation.

Including the 5th term has a very specific value. It makes the error distribution look quite Gaussian. By our engineer’s eyes (or hearts), it should be very desirable. Outliers can be very costly in an engineering context, leading to delays and frustration. Just to close the subject, by outlier, I mean an event whose frequency is not predicted by Gaussian statistics, such as an error of 6 standard deviations occurring once in a sample of a thousand. Such events should require two billion of samples.

Commentaires sur le texte

Il y a des choses dans la scénario qui méritent qu'on s'y attarde un peu. L'une des plus troublantes est que lorsque nous interpolons et rééchantillonnons les vecteurs singuliers gauche et droit (colonnes des matrices U et V), ils ne sont plus orthogonaux. Cela signifie que nous ne pouvons pas utiliser les vecteurs interpolés pour projeter d'autres distributions sur cet espace. Le manque d'orthogonalité ne semble cependant pas ajouter de biais à la reconstruction.

En outre, les interpolations utilisent le nombre de valeurs singulières extraites comme ordre du polynôme utilisé pour l'interpolation. Il est facile d'argumenter en faveur de ce choix pour un modèle comme celui-ci. Ce choix, cependant, serait l'un des premiers que j'essaierais de confirmer (ou de rejeter) lorsque je serais confronté à des données réelles.

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.

Merci de votre lecture.