Distortion Mapping – Un HowTo per l'SVD

In questo articolo tecnico vorrei discutere di come (e probabilmente come non farlo) sviluppare una correzione delle distorsioni ottiche. problema ottico comune, e uno che in realtà può essere compensato utilizzando un algoritmo perché, anche se l'immagine è distorta, può comunque essere priva di aberrazioni frontali. Credo che il titolo dica già in volgarità che si tratterà di un seguito della mia storia d'amore con la decomposizione del valore singolare.

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.

Assunti sul problema da risolvere

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.

Il problema della correzione del problema della distorsione si inserisce bene nel quadro teorico della SVD perché possiamo raccogliere i dati per adattarli al formalismo della matrice. Pertanto, supporò che raccoglieremo le informazioni sulla distorsione su una griglia che copra uniformemente un'immagine rettangolare.

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 qui. (Siete invitati a inviare una e-mail a Contatto 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

Il codice sopra è solo un'espansione di MacLaurin di ordine 5 con coefficienti casuali. L'espansione definisce i spostamenti lungo l'asse delle coordinate su un griglia che analizziamo utilizzando la decomposizione dei valori singolari. Il risultato di un'analisi dei valori singolari è una decomposizione della matrice in tre termini,

The matrix “S” is diagonal and (conventionally) sorted in descending order. The columns of the U e V Le matrici sono ortogonali e (di solito) normalizzate,

Prima di approfondire i dettagli dell’analisi della distorsione del nostro modello, vorrei collegare la SVD a qualcosa che la maggior parte dei lettori potrebbe già conoscere: il metodo della separazione delle variabili. Quando vogliamo risolvere analiticamente una equazione differenziale parziale, assumiamo che la soluzione possa essere scritta come:,

Questo ci porta alla soluzione che scriviamo come:,

Questo è esattamente lo scopo della decomposizione del valore singolare: offre un modo per scrivere la nostra 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 e poteri y multiplicata per un coefficiente (random) specifico. L'analisi SVD ha ridotto il nostro insieme di 20 funzioni input a due set (U e V) di 6 funzioni. Questo si può vedere eseguendo lo script (una versione completa dello script può essere scaricata). qui)

I risultati sopra riportati sono stati ottenuti senza aggiungere rumore, per distinguere chiaramente ciò che è una soluzione e ciò che è (in numeri) uguale a zero.

Analizzando la distorsione con il rumore

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,

Con l'aggiunta del rumore, la struttura dei valori singolari appare diversa. Mentre nel primo caso i valori singolari diminuiscono di 12 ordini di grandezza dopo i primi 6, ora si ha quello che sembra una sovrapposizione tra due linee rette. In realtà questo è abbastanza comune quando si aggiunge rumore bianco. Per il nostro modello di distorsione, il rumore bianco è probabilmente una buona ipotesi perché le misure dovrebbero essere indipendenti.

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.

Commenti sullo script

Ci sono alcune cose nel sceneggiatura Questo aspetto potrebbe valere la pena di essere approfondito. Ciò che mi sembra più inquietante è che quando interpoliamo e ridimensioniamo i vettori singolari sinistri e destri (le colonne delle matrici U e V), questi non sono più ortogonali. Ciò significa che non possiamo utilizzare i vettori interpolati per proiettare altre distribuzioni su questo spazio. La mancanza di ortogonalità non sembra tuttavia comportare alcun pregiudizio nella ricostruzione.

Inoltre, le interpolazioni utilizzano il numero dei valori singolari estratti come ordine del polinomio utilizzato per l’interpolazione. È facile sostenere questa scelta per un modello del genere. Tuttavia, questa scelta sarebbe una delle prime che vorrei provare a confermare (o a rifiutare) quando mi trovo di fronte a dati reali.

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.

Grazie per aver letto.