Antwort auf Frage 2 aus Ihrem Qualifikationstest.
Von : Menno Rubingh
Datum: 21.08.2011
--------------------------------------------------------------------------
NB: Dieser Text is ISO-8859-1 plain text, formatiert für viewing mit einem
MONOSPACE font.
--------------------------------------------------------------------------
Die Literatur die ich gesehen habe seit Fr 12.08. ist:
Trucco & Verri, 1998,
_Introductory Techniques for 3-D Computer Vision_
Zhang, 2000,
das Paper aus IEEE Transactions on Pattern Analysis and Machine
Intelligence, Vol 22, No 11 (p. 1330 - 1334).
Cyganek & Siebert, 2009,
_An Introduction to 3D Computer Vision Techniques and Algorithms_
ISBN 978-0-470-01704-3
Wöhler, 2009,
_3D Computer Vision: Efficient Methods and Applications_
ISBN 978-3-642-01731-5
Kurze Schilderung meiner Vorkenntnisse vor dem 12.08.:
Vor dem 12.08. war das Thema "camera calibration" mir noch völlig
unbekannt. In meinem Job 2008-2010 bezüglich 3D CAD, für DAKO GmbH in
Jena, hatte ich jedoch schon gelernt über einfache projective geometry
(pinhole-Kamera), homogene Koordinaten, und das mapping 2D conic <-->
3D Parabol. (Epipolar lines und den "absolute conic" waren mir noch
unbekannt.)
Fürs Auffrischen meiner Kenntnisse Über singular value decomposition
und Levenberg-Marquard habe ich nachschlagen müssen in Strang _Linear
Algebra and its Applications_; Stoer&Bulirsch _Introduction to Numerical
Analysis_; und Press et.al. _Numerical Recipes_.
--------------------------------------------------------------------------
Inhalt
* KAMERAMODELL
* KAMERAMODELL IN HOMOGENEN KOORDINATEN
* KALIBRIERUNG
- Kalibrierung aus N Korrespondenzen, ohne radial distortion
- Verfeinerung der Kamera-parameters
- Planarer calibration rig
===========================================================================
KAMERAMODELL
===========================================================================
Ein "Kameramodell" beschreibt wie die 3D Position eines "world
Punktes" (in der fotografierten Szene) sich abbildet zu (maps to) den
pixel-Koordinaten in dem digitalen Bild das aus der Kamera ausgelesen
wird.
Der wesentlichste Schritt dieses mapping (Abbildung) ist das mapping von
3D Punkt zur der 2D Flaeche worin der Bildsensor liegt. Das modellieren
der Lichtfall (deflection) durch eine optischen Linse (lens) ist komplex;
man benutzt das angenäherte (Approximate) Modell eines "pinhole camera":
xP = - f x1 / z1 (Gleichung F1a)
yP = - f y1 / z1 (Gleichung F1b)
Parameter f = "focal length".
x1,y1,z1 = 3D Koordinaten des scene point.
xP,yP = 2D Koordinaten in Bild-flaeche.
Das x1,y1,z1 Koordinatensystem (= "Kamera-Koordinatensystem"):
Z-Achse perpendicular zur Bild-Flaeche, und
Ursprung in einem Z-Abstand f von der Bildfläche.
Das xP,yP Koordinatensystem:
Ursprung = wo die Z-Achse des x1,y1,z1 Koordinatensystems die
Bildfläche schneidet (= "principal point");
X- und Y-Achsen gleich wie im x1,y1,z1-Koordinatensystem.
Die Minus-Zeichen weil das Bild "umgekehrt" steht.
[ Jeder 3D Punkt ( w x1, w y1, w z1 ) auf einer Linie durch (x1,y1,z1)
und den Ursprung des Kamera-Koordinatensystems wird abgebildet auf den
gleichen 2D Punkt (xP,yP). ]
Wenn die Abweichungen relativ zur pinhole-Annäherung, verursacht durch
das Linsensystem, zu groß sind, dann ist die nächste Verfeinerung dass
man die folgende Korrektur ("radial distortion") benutzt:
alfa xPN = xPV
alfa yPN = yPV
wo
alfa = 1 + k1 r^2 + k2 r^4
mit
r^2 = xPN^2 + yPN^2
k1 = Parameter
k2 = Parameter (= 0 in erster Annäherung)
xPV,yPV = 2D Koordinaten in Bild-flaeche, VOR Korrektur
xPN,yPN = 2D Koordinaten in Bild-flaeche, NACH Korrektur.
[ NB: Hier scheint es mir dass man diese Korrektur ausführen sollte auf
xP,xP -- und nicht auf xS, yS wie behauptet in (Cyganek & Siebert,
2009). ]
Dann muss man die (eventuell korrigierten) xP,yP umrechnen zu den
Pixel-Koordinaten xS,yS des Bildsensors (= des digitalen Bildes):
xP = Hx ( xS - Ox )
yP = Hy ( yS - Oy )
Parameter (Ox,Oy) = Position des principal point in Pixel-Koordinaten.
Parameter Hx,Hy = Breite und Höhe eines Pixels.
Das xS,yS Koordinatensystem:
Ursprung in einer Ecke des Bildes.
xS und yS müssen hier nicht unbedingt integer Werte haben.
NB: Jedes der Koordinatensysteme {x1,y1,z1}, {xP,yP}, {xS,yS} hat ihre
X-Achse und Y-Achse parallel zu den Reihen und Spalten der Pixel.
Die bis soweit genannten Parameter nennt man "intrinsic parameters",
weil sie mit der Kamera selbst zu tun haben, und nicht mit der Position
u/o Orientierung der Kamera in der fotografierten 3D-Welt. Die intrinsic
parameters sind also:
f = focal length
k1, k2 = radial distortion.
Ox, Oy = Position des principal point in Pixel-Koordinaten.
Hx, Hy = Breite und Höhe eines Pixels.
Wenn man die radial distortion Korrektur nicht ausführt -- ODER auch
wenn man die Korrektur (m.E. im Grunde fälschlicherweise) durchführt
auf xS,yS anstatt auf xP,yP -- , dann ergibt das zusammennehmen der
obigen Formeln für xS,yS (vor eventuellem Korrektur):
xS = - (f/Hx) (x1/z1) + Ox
yS = - (f/Hy) (y1/z1) + Oy
worin man sieht dass dass f und Hx immer zugleich auftreten als Faktor
(f/Hx); und genauso für f und Hy. Das bedeutet dass die drei Parameter
{ f, Hx, Hy } in diesem Fall nur zwei degrees of freedom haben.
In diesem Fall würde ich anstatt { f, Hx, Hy } die zwei folgenden
Parameter benutzen:
fx = f / Hx = ratio ( focal length / pixel size ) für X
fy = f / Hy = ratio ( focal length / pixel size ) für Y.
Jedoch wenn man die Korrektur macht auf xP,yP, dann würde das
bedeuten dass f, Hx, Hy dann im Kameramodell NICHT immer nur als
(f/Hx) und (f/Hy) vorkommen. Die Parameter { f, Hx, Hy } würden
dann 3 degrees of freedom haben, was sich jedoch nur in der radial
distortion äussert: was m.E. bedeutet dass die Kalibrier-
Berechnungen dann nur funktionieren bei einer Kamera die
tatsächlich eine signifikante radial distortion hat, und auch nur
dann wenn sich genug Korrespondenz-Punkte in einer stark
distortierten Zone im Bild befinden.
Am besten scheint es mir, dass man beim Kalibrieren aus den Daten
detektiert ob radial distortion signifikant ist oder nicht. Wenn
signifikant, dann xP,yP auf radial distortion korrigieren und die
intrinsic Parameter { f, k1, k2, Ox, Oy, Hx, Hy } benutzen (wobei man
k2 = 0 nehmen sollte wenn k2 nicht signifikant ist).
Wenn distortion nicht signifikant ist, dann die intrinsic Parameter
{ fx, fy, Ox, Oy } benutzen, und die radial distortion Korrektur im
Modell disablen (k1 = k2 = 0).
Die rotation und translation des Kamera-Koordinatensystems x1,y1,z1
relativ zu einem anderen Koordinatensystem W in der 3D Welt nennt man
die "extrinsic parameters":
p1 = R ( pW - t )
wo
p1 = (x1,y1,z1)
= Koordinaten des 3D Punktes im Kamera-Koordinatensystem
pW = (xW,yW,zW)
= Koordinaten des selben 3D Punktes im W-Koordinatensystem
t = translation vector = (tX,tY,tZ).
R = 3x3 Rotationsmatrix (orthonormal).
Die Rotationsmatrix R hat nur drei unabhängige Parameter (die
drei Rotationswinkel um die Koordinat-Achsen; oder Richtung
der Rotationsachse plus Rotationswinkel). Zusammen mit dem
Translationsvektor t ergibt das insgesamt 6 unabhängige extrinsic
parameters.
===========================================================================
KAMERAMODELL IN HOMOGENEN KOORDINATEN
===========================================================================
Im oben dargelegten Kameramodell ist die Beziehung zwischen 3D Welt-Punkt
(xW,yW,zW) und 2D Bild-Punkt (xS,yS) nicht linear, einerseits wegen
der fundamentalen Gleichungen F1a,F1b für die Perspectiv-Abbildung,
und andererseits wegen der radial distortion.
Ein sehr häufig benutzter Trick, das Kameramodell trotzdem in linearen
Gleichungen zu fassen, besteht darin dass man
- radial distortion ausser Acht lässt [note 1], und
- sogenannte "homogene Koordinaten" einsetzt für das 2D Bild-Punkt.
Da ich die über diesen Trick gewonnenen linearen Gleichungen benutze
im Kapitel unten über Kalibrierung, fasse ich hier kurz zusammen wie
der Trick funktioniert.
Bei ignorieren von radial distortion ist die Beziehung
zwischen Bild-Punkt (xS,yS) und einen 3D-Punkt (x1,y1,z1) im
Kamera-Koordinatensystem die folgende:
xS = - fx (x1/z1) + Ox
yS = - fy (y1/z1) + Oy
Durch einführen von neuen Variabelen xS~, yS~, zS~ (genannt "homogene
Koordinaten", und angezeigt mit einer '~') kann man dies úmschreiben als
xS = xS~ / zS~ (Gleichung M1a)
yS = yS~ / zS~ (Gleichung M1b)
wo
xS~ = -fx x1 + Ox z1
yS~ = -fy y1 + Oy z1
zS~ = z1
Die Beziehung zwischen den xS~, yS~, zS~ (woraus man immer direkt die
xS,yS finden kann) und p1 = (x1,y1,z1) ist also ein System von linearen
Gleichungen:
+- -+ +- -+ +- -+
| xS~ | | -fx 0 Ox | | x1 |
| yS~ | = | 0 -fy Oy | | y1 | (Gleichung M2)
| zS~ | | 0 0 1 | | z1 |
+- -+ +- -+ +- -+
\_____________/
Mint
Die Matrix mit den intrinsic parameters in Gleichung M2 nennt man: Mint.
Auch das das mapping zwischen p1 = (x1,y1,z1) im Kamera-Koordinatensystem
und pW = (xW,yW,zW) im W-Koordinatensystem kann man schreiben als System
von linearen Gleichungen. Aus
p1 = R ( pW - t )
macht man zunächst
p1 = R pW - tR
wo tR = R t = der rotierte translation vector. (Translation t und
dann Rotation R ist gleich zu: Die gleiche Rotation R und dann eine
Translation um den rotierten translation vector tR.) Letztere Gleichung
ausgeschrieben in Koordinaten ist:
+- -+ +- -+ +- -+
| r11 r12 r13 | | xW | | tRx |
p1 = | r21 r22 r23 | | yW | - | tRy |
| r31 r32 r33 | | zW | | tRz |
+- -+ +- -+ +- -+
Das ist das gleiche wie:
+- -+ +- -+
| r11 r12 r13 -tRx | | xW |
p1 = | r21 r22 r23 -tRy | | yW |
| r31 r32 r33 -tRz | | zW | (Gleichung M3)
+- -+ | 1 |
+- -+
\____________________/
Mext
wo der Vektor (xW,yW,zW) um einen Koordinat erweitert ist, der immer
den Wert 1 hat. Gleichung M3 ist das gesuchte lineare Gleichungssystem
für die Beziehung zwschen p1 und pW. Die Matrix mit den extrinsic
parameters in Gleichung M3 nennt man: Mext.
Insgesamt ist die Beziehung zwischen pW = (xW,yW,zW) und pS also
xS = xS~ / zS~ (Gleichung M1a)
yS = yS~ / zS~ (Gleichung M1b)
+- -+ +- -+
| xS~ | | xW |
| yS~ | = M | yW | (Gleichung M4)
| zS~ | | zW |
+- -+ | 1 |
+- -+
wo M = Mext Mint.
Ausgeschrieben:
+- -+ +- -+ +- -+ +- -+
| xS~ | | -fx 0 Ox | | r11 r12 r13 -tRx | | xW |
| yS~ | = | 0 -fy Oy | | r21 r22 r23 -tRy | | yW |
| rS~ | | 0 0 1 | | r31 r32 r33 -tRz | | zW |
+- -+ +- -+ +- -+ | 1 |
+- -+
\____________/ \____________________/
Mint Mext
(Gleichung M5)
===========================================================================
KALIBRIERUNG
===========================================================================
"Kalibrierung" bedeutet dass man die Werte der Parameter eines
Kamera-Modells feststellt aus Messungen (measurements).
Die einfachste Art der Kalibrierung ist dass man eine 3D Szene
fotografiert die eine Anzahl 3D Punkte enthält deren W-Koordinate
(xW,yW,zW) mann kennt. Diese 3D-Punkte sollten erkennbar sein im
aufgenommenen Bild. Man findet dann die Korrespondenzen welcher Punkt
pS im Bild zu welchem Welt-Punkt pW gehört.
Jede Korrespondenz ist ein Paar { pW, pS }, wo pW = (xW,yW,zW) und pS =
(xS,yS) bekannt sind. Für jede Korrespondenz schreibt man dann die
Formeln des Kamera-Modells auf. In diesem System von Gleichungen sind
alle xW,yW,zW und alle xS,yS bekannt, und alle Parameter (intrinsic und
extrinsic parameters) unbekannt. Wenn man genug Korrespondenzen hat,
kann man hieraus die Werte der Parameter finden. (Mehr hierüber unten.)
Ergebnis ist der Satz der Parameter-Werte.
[ Nebenbei: Die extrinsic parameters (R, t) sind natürlich nur gültig
solange man die Kamera nicht verschiebt/dreht im W-Koordinatensystem
(der 3D-Kalibrierpunkte). D.h. die extrinsic parameters
muss man neu bestimmen (neu vermessen) sobald die Kamera
umpositioniert/gedreht wird.
Verschieben/drehen der Kamera sollte aber die intrinsic
parameters (für eine gegebene Kamera) nicht beeinflussen, solange
Faktoren konstant bleiben die einen Einfluss haben könnten auf
die Parameter-Werte (z.B. Temperatur). ]
Die Information die man pro Korrespondenz aus den Messdaten gewinnt
ist die folgdende. Das Kamera-Modell gibt pro Korrespondenz zwei
Gleichungen:
xP = - f (x1/z1)
yP = - f (y1/z1)
[ wo (xP,yP) von (xS,yS) abhängt, und (x1,y1,z1) von (xW,yW,zW)
abhängt; aber Substitution dieser Beziehungen in obige Gleichungen
ändert nicht die Anzahl dieser obigen Gleichungen ].
Wenn man N Korrespondenzen einbezieht, sind das insgesamt 2N Gleichungen,
in denen die N 3D-Kalibrierpunkte (xW,yW,zW) und deren N 2D-Abbildungen
(xS,yS) bekannt sind, und in denen die Kamera-parameters unbekannt sind.
Theoretisch sollte man daher aus N Korrespondenzen die Werte von
_maximal_ 2N Kamera-parameters feststellen können.
Das Kamera-Modell ohne radial distortion hat die 4 intrinsic parameters
{ fx, fy, Ox, Oy }, also mit den 6 extrinsic parameters insgesamt 10
parameters. Dafür würde man mindestens 10 Gleichungen brauchen, d.h.
mindestens 5 Korrespondenzen.
Mit radial distortion hat das Kameramodell 7 intrinsic parameters
{ f, k1, k2, Ox, Oy, Hx, Hy }, also also mit den 6 extrinsic parameters
insgesamt 13 parameters; wofür mindestens 13 Gleichungen benötigt sind,
d.h. mindestens 7 Korrespondenzen.
Wie bei jedem Messverfahren, erreicht man eine bessere Messung wenn
man mehr Messpunkte (d.h. mehr Korrespondenzen) berücksichtigt.
Eine andere Art von Kalibrierung ist die sogenannte "self-calibration".
Dabei ist es nicht nötig dass die Koordinaten bekannt sind der 3D
Punkte die man verwendet (und im Bild erkennt). Jedoch man braucht
dabei von dieser 3D Szene (die die 3D Punkte enthält) mehrere Aufnahmen,
aus unterschiedlichen Kamerapositionen. (Die Position der verwendeten
3D Punkte im W-Koordinatensystem sollte sich zwischen den Aufnahmen
nicht ändern.)
Bei self-calibration besteht die Information die man ausnutzt um
die Kamera-parameters (und die unbekannte 3D position der erkannten
3D scene points) festzustellen aus: Die Korrespondenzen zwischen (a)
der 2D-Abbildung in Kameraposition 1 der erkannten K 3D scene points,
und (b) dessen 2D-Abbildung in Kameraposition 2.
Ich behandle in diesem Text nur Kalibrierung auf grund von 3D
Kalibrierpunkten mit BEKANNTEN Koordinaten im W-Koordinatensystem,
und nicht self-calibration.
In den ersten zwei Abschnitten unten beschreibe und erkläre ich
eine genaue Vorgehensweise fürs ableiten (feststellen) der Werte der
Modellparameter aus N Korrespondenzen zwischen bekannten 3D Punkten
und Ihre 2D-Abbildungspunkten, in einer eizigen Aufnahme:
- Abschnitt "Kalibrierung aus N Korrespondenzen, ohne radial distortion"
Hier beschreibe ich wie man aus den Messdaten _analytisch_ eine
Schätzung der Kamera-parameters finden kann. Der Algorithmus
den ich hier beschreibe und analysiere ist der aus Trucco&Verri,
Abschnitt 6.3.1 (Seite 132-133). Dieser Algorithmus benutzt das
Kameramodell in homogenen Koordinaten (Gleichung M5 oben), also
ohne radial distortion (k1= k2 = 0).
- Abschnitt "Verfeinerung der Kamera-parameters"
Hier erkläre ich wie man die analytisch gefundene Schätzung der
Kamera-parameters verfeinern (genauer machen) kann, und auch
erweitern kann um das finden der radial distortion parameters.
Dies geschieht aus den gleichen Messdaten, und über eine numerische
(iterative) Prozedur in der das volle, nicht-lineare, Kameramodell
benutzt wird.
Die iterative Verfeinerungs-Prozedur braucht einen initial guess.
Als initial guess wird benutzt: k1 = k2 = 0, plus die angenäherten
Kamera-parameters gefunden mit der analytischen Methode im ersten
Abschnitt.
Insgesamt ist die Vorgehensweise also dass man die Schritte aus den
zwei obengenannten Abschnitten beide ausführt, nach einander. Also
zunächst analytisch einen initial guess finden, danach diesen verfeinern.
Im ersten Abschnitt erkläre ich anhand der dort beschriebenen
Kalibriervorgang auch, wie viele 3D Kalibrierpunkte benötigt sind,
und welche Bedingungen die 3D Kalibrierpunkte weited erfüllen müssen.
Eine wichtigste Bedingung die dort "entdeckt" (erklärt) wird ist
dass die 3D Kalibrierpunkte nicht alle in der gleiche Fläche (plane)
liegen sollten. Diese Bedingungen gelten _auch_ für den zweiten
Abschnitt ("Verfeinerung").
Es gibt jedoch Verfahren zur Kamerakalibration die calibration rigs
(Kalibriermuster) benutzen die PLANAR sind, wie das Verfahren von
Zhang 2000. Diese Verfahren behandelt der letzte Abschnitt:
- Abschnitt "Planarer calibration rig"
Hier analysiere ich welche Kalibrierungen möglich sind mit einem
planaren Kalibriermuster.
Hauptteil des Abschitts ist meine Analyse von dem Verfahren von
Zhang 2000 für das Feststellen der (intrinsic) Kamera-parameters
mittels eines planaren Kalibriermusters, wovon _mehrere_ Aufnahmen
gemacht werden. Ich erkläre dabei welche Bedingungen diese
Aufnahmen erfüllen müssen (nämlich: Man sollte den calibration rig
zwischen Aufnahmen verschieben und/oder drehen relativ zur Kamera).
Alle Kalibrier-Methoden die ich in diesem Text beschreibe beruhen darauf,
dass man in einem Vorgang zugleich alle Kamera-Parameter bestimmt.
Ich habe mich auf diese Art der Kalibration konzentriert weil auch
Zhang 2000 so vorgeht.
Jedoch man sollte m.E. nicht vergessen dass es auch Methoden gibt,
die einzelne intrinsic parameters direkt festlegen, d.h. über Verfahren
die nicht zugleich auch andere intrinsic parameters vermessen. In der
gesehenen Literatur fand die folgende erwähnt:
- Ox,Oy :
Trucco&Verri, Abschnitt 6.2.3 (p.131-132) bestimmen die Pixel-
Koordinaten Ox,Oy des principal point aus den vanishing points
einer Aufnahme eines Kalibriermusters das mehrere unterschiedliche
Mengen gerader Linien in 3D festlegt.
- Radial distortion:
Wöhler p. 18 schreibt laconisch: "Radial and tangential distortions
can be directly inferred from deviations from straightness in the
image" (aber beschreibt dort selbst keine detaillierte Methode).
----------------------------------------------------------
Kalibrierung aus N Korrespondenzen, ohne radial distortion
----------------------------------------------------------
>>> Abfassen der N Korrespondenzen als Matrixgleichung A m = 0
Wie erklärt im Kapitel "KAMERAMODELL IN HOMOGENEN KOORDINATEN" oben, kann
man, wenn man radial distortion ausser Acht lässt, homogene Koordinaten
einsetzen, und damit das Kamera-Modell schreiben als lineare Gleichungen:
+- -+ +- -+
| xS~ | | xW |
| yS~ | = M | yW | (Gleichung M4)
| zS~ | | zW |
+- -+ | 1 |
+- -+
wobei
xS = xS~ / zS~ (Gleichung M1a)
yS = yS~ / zS~ (Gleichung M1b)
Die homogene Koordinaten sind mit einer ~ versehen. Die Gleichungen M1a,
M1b definieren das mapping zwischen den homogenen Koordinaten xS~,yS~,zS~
und die gemessenen Bild-Koordinaten xS,yS.
M ist ein 3x4 Matrix, in der die Elemente m_ij lineare Kombinationen
der Parameter sind (siehe Gleichung M5 weiter oben). Die Gleichungen
M1a, M1b zeigen dass man den right hand side in Gleichung M4 mit einem
willkürlichen scalar factor multiplizieren kann ohne dass xS,yS sich
ändern; d.h. Matrix M ist definiert up to an arbitrary multiplicative
factor.
Obiges System vom Gleichungen gibt es für jede Korrespondenz.
Die bekannten Messwerte sind xS,yS und xW,yW,zW (natürlich andere
Werte für jede Korrespondenz). Unbekannt (und zu bestimmen) sind die
Kamera-Parameter in M (die sind die gleichen für jede Korrespondenz).
Trucco&Verri 1998, p. 133, gehen jetzt so vor dass zunächst die Elemente
m_ij der Matrix M als die unknowns behandelt werden. (Hieraus werden
dann später die Kamera-Parameter berechnet.) Sie schreiben die
Gleichungen M1a,M1b so dass sie linear werden in den unknowns m_ij:
xS~ - xS zS~ = 0
yS~ - yS zS~ = 0
Substitution hierin der xS~, yS~, zS~ durch die right hand sides aus
Gleichung M4 ergibt zwei Gleichungen die linear sind in den unknowns
m_ij. Für N Korrespondenzen ergibt das N dieser Gleichungspaare (alle
linear in den unknowns m_ij). Jetzt ist erreicht die Matrixgleichung
(6.17) aus Trucco&Verri 1998, Seite 133:
A m = 0 (Gleichung K1)
wo
m = vector der 12 Elemente der Matrix M
= ( m_11, m_12, m_13, m_14, m_21, m_22, ... , m_34 )
und wo A der folgende 2N x 12 Matrix ist:
+- -+
| xW1 yW1 zW1 1 0 0 0 0 -xS1*xW1 -xS1*yW1 -xS1*zW1 -xS1 |
| 0 0 0 0 xW1 yW1 zW1 1 -yS1*xW1 -yS1*yW1 -yS1*zW1 -yS1 |
A = | ... ... ... |
| xWN yWN zWN 1 0 0 0 0 -xSN*xWN -xSN*yWN -xSN*zWN -xSN |
| 0 0 0 0 xWN yWN zWN 1 -ySN*xWN -ySN*yWN -ySN*zWN -ySN |
+- -+
(Gleichung K2)
>>> Lösen von A m = 0,
und Bedingungen für die 3D Kalibrationspunkte
Aus Gleichung K1 (A m = 0) sehen wir noch einmal dass die m_ij nur
definiert sind up to an arbitrary multiplicative factor. Mathematisch
gesehen sagt Gleichung K1 dass der gesuchte Vektor m ein beliebiger
Vektor ist im nullspace von Matrix A.
"Theoretisch" (ohne Messfehlern, ohne Diskretisierungsfehlern, ohne
Abrundungsfehlern, usw., und mit 100% Übereinstimmung zwischen Modell und
Wirklichkeit), und wenn N groß genug ist, kann die 12-spaltige matrix
A maximal den rank 11 haben (wobei die nullspace die Dimension 1 hat).
Rank 11 reicht aus dazu, die 11 Freiheitsgrade in m festzulegen.
Damit A den rank 11 erreicht sollte sie natürlich mindestens 11 Zeilen
(rows) enthalten; da jede Korrespondenz 2 rows beiträgt brauchen wir
also mindestens 6 Korrespondenzen. [note 2]
Darüber hinaus gibt es noch eine Bedingung damit A den rank 11 erreicht,
nämlich: Es sollte nicht so sein dass alle 3D Kalibrierungspunkte pW_i =
(xWi,yWi,zWi), i = 1 ... N, in der gleichen Fläche (plane) liegen.
Betrachten wir z.B. die Fläche zW = 0. [Das ist allgemein genug, da
die Orientierung des W-Koordinatensystems (= Rotation R) jede sein kann.]
Bei Substitution von zWi = 0 für alle i in Matrix A (Gleichung K2)
entstehen 3 null-wertige Spalten. Der rank von Matrix A verringert
sich hierdurch um 3, also von 11 nach 9. (Der nullspace von A ist jetzt
nicht mehr 1-dimensional, sondern 4-dimensional.)
Diese 9 freiheitsgrade reichen nicht aus dazu, das Kameramodell zu
kalibrieren dass neben den 6 extrinsic parameters die 4 intrinsic
parameters { fx, fy, Ox, Oy } benutzt.
Mehrere der gelesenen Literaturstücke schreiben dass Ox = Oy = 0
für viele Kameras eine gute erste Annäherung ist. Das verringert die
Anzahl intrinsic parameters zu 2, nämlich { fx, fy }, also mit den
extrinsic parameters insgesamt 8 Kamera-parameters. Dieses Kameramodell
würde man kalibrieren können aus einer einzigen Aufname eines planaren
Kalibriermusters; siehe dazu Abschnitt "Planarer calibration rig" unten.
In der Wirklichkeit gibt es immer Messfehler, Diskretisierungsfehler, und
Abrundungsfehler, und auch benutzen wir ein vereinfachtes (angenähertes)
Kameramodell. Schon bei 6 Korrespondenzen wird die matrix A daher
in Wirklichkeit nicht den rank 11 haben, sondern den vollen rank 12.
Der nullspace ist jetzt verschwunden und das Gleichungssystem K1 hat
keine Lösung mehr (außer die für uns nutzlose triviale Lösung m = 0).
Wir können trotzdem weitermachen und einen _optimalen_ Vektor m
bestimmen, dadurch dass wir die Fehler (Abweichungen wegen Messfehlern
usw.) modellieren als
A m = e, e = Fehlervektor
und dann der Vektor m suchen der die Länge || e || der Fehlervektors
minimiert. (Dies ist least squares Optimierung: Es wird minimiert der
sum of squares der Elemente des Fehlervektors.) Diese Vorgehensweise
erlaubt es auch, die Anzahl N der Korrespondenzen zu vergrößern, so
dass matrix A höher wird als 12 Zeilen. Minimieren von e bedeutet dann
dass der gefundene m ein "fit" ist des Modells auf alle N Messpunkte.
Je mehr Messpunkte einbezogen werden, desto genauer sollte der gefundene
m sein (die Annahme dabei ist jedoch dass die Abweichungen über alle
Messpunkte statistisch gleich verteilt sind).
Der m der || e || = || A m || minimiert findet man m.E. am besten aus
Singular Value Decomposition (SVD) A = U D V^T der Matrix A.
|| A m || ist minimal wenn m = v1, wo v1 die Spalte aus V ist die gehört
zur kleinsten singular value s1 (aus D).
Je kleiner die Abweichungen e in den Messdaten sind (und je mehr
Modell mit Wirklichkeit übereinstimmt), desto mehr nähert sich der
kleinste singular value zu 0 an. Im theoretischen Idealfall e = 0 ist
der kleinste singular value s1 gleich null, und liegt m = v1 in der
nullspace und ist ein exacter fit (mit || e || = 0).
Der kleinste singular value s1 ist also m.E. ein Maß für die "Güte"
des least squares fit: Je kleiner desto besser (0 = ideal); und ein
großer s1 bedeutet dass Messdaten und Modell nur schlecht übereinstimmen.
Es sollte weiter so sein dass A nur _einen_ singular value hat der nah
bei 0 liegt; die andere singular values sollten "deutlich" größer als
0 sein. Wenn es (trotz N >= 6) mehrere near-zero singular values gibt,
so dass Richtung m die || A m || minimiert schlecht festgelegt ist, dann
bedeutet das, dass die 3D Kalibrierungspunkte nahezu in der gleichen
Fläche liegen (vergleiche die Analyse des Idealfalls oben).
M.E. sollte die Kalibrierungs-Berechnung die in den letzten Absätzen
genannten Bedingungen überprüfen, und es dem Anwender melden wenn die
Bedingungen schlecht erfüllt sind. (Man könnte z.B. dem Anwender die
folgenden Werte zurückgeben: Der kleinste singular value s1; die ratio
r1N = s1/sN des kleinsten und grössten singular value; und die ratio
r12 = s1/s2 der zwei kleinsten singular values.)
Zhang 2000 schreibt zu recht noch dass ein Matrix wie unsere A
(Gleichung K2) numerisch schlecht konditioniert ist: Die Spalten 1,2,3
und 5,6,7 enthalten W-Koordinaten, die Spalten 4 und 8 enthalten units,
die letzte Spalte enthält Bild-Koordinaten, und die Spalten 9,10,11
enthalten sogar Produkte von W-Koordinaten und Bild-Koordinaten --
d.h. die Spaltenvektoren von A sind widely out of scale. Ich würde
hier vor SVD jede Spalte von A normalisieren auf die gleiche Länge
wie die Spalten 4 und 8. Die benutzte Normierungsfaktoren für jede
Spalte speichern, und nachher die elemente m_ij des gefundenen Vektors
m zurückskalieren um den vorher für die entsprechenden Spalte aus A
benutzten Faktor.
Nachdem die m_ij erfolgreich festgestellt sind, berechnet man daraus
die intrinsic parameters, die elemente des rotation matrix, und den
translation vector tR, über analytische Formeln die sich ohne Probleme
(siehe Trucco&Verri, p. 134 - 135) ableiten lassen aus M = Mint Mext,
mit den Mint und Mext aus Gleichung M5, und wo M die jetzt bekannten
m_ij enthält (siehe Trucco&Verri, p. 134 - 135).
Die Elemente so erhaltene rotation matrix R sind aus Messungen
abgeleitet, und diese R is daher wahrscheinlich nicht orthogonal,
und sollte korrigiert werden, wie folgt: SVD von R berechnen (R =
U D V^T). Die korrigierte rotation matrix ist dann U I V^T, mit den
singular values auf 1 gesetzt.
----------------------------------
Verfeinerung der Kamera-parameters
----------------------------------
Wöhler 2009 (u.a. um Seite 22) kritisiert, dass der error Vektor
e = A m, der minimiert wird im obigen Abschnitt, eine arbiträre Größe
ist, und keine physische Bedeutung hat.
Und er schreibt dass die Parameter-Werte die man berechnet aus den so
erhaltenen m_ij nur taugen als initial guess in einer Optimierung in
der man stattdessen den "reprojection error" E minimiert:
E = SUM_i [ d_i ^ 2 ] (Gleichung K3)
mit
d_i = abstand zwischen Punkt pS_i und Punkt Pr( pW_i )
wo
pS_i = Das 2D Bild-Punkt der i-ten Korrespondenz
(in pixel coordinates)
pW_i = Das 3D Welt-Punkt der i-ten Korrespondenz
Pr( pW_i ) = Aus Kamera-Modell Berechnete Projektion von
pW_i nach 2D Bild-Punkt
(in pixel coordinates).
In der Berechnung der Pr( pW_i ) aus den pW_i fliessen hinein die
Werte der Kamera-parameters. Minimieren von E bedeutet daher dass
man die Kamera parameters so sucht dass die Pr( pW_i ) nah zu ihren
korrespondierenden pS_i liegen.
(Auch von Zhang 2000 wird die gleiche error measure E minimiert wenn
er seine Verfeinerung ausführt.)
Pr( pW ) ist eine nicht-lineare Funktion der Koordinaten von pW =
(xW,yW,zW). Minimieren von E (als Funktion der Parameter) bedeutet daher
nicht-lineare Minimierung, d.h. Levenberg-Marquard [beshrieben in Press
et.al., _Numerical Recipes_]. Jeder Algorithmus für nicht-lineare
Minimierung braucht eine initial guess (in unserem Fall: für die
Kamera-Parameter) -- d.h. die nicht-lineare Minimierung *verfeinert*
die initial guess.
Das alles passt neatly zusammen: Zunächst findet man aus linearen
Gleichungen (mit radial distortion ignoriert) über SVD und über
minimieren einer "komischen" error measure eine initial guess für die
Kamera-parameters. Diese initial guess verfeinert man dann mit ein paar
Levenberg-Marquard Schritten zur minimierung der richtigen error measure
E, wobei man auch gleich das Modell erweitern kann um radial distortion
(initial guess: k1 = k2 = 0).
------------------------
Planarer calibration rig
------------------------
Der Algorithmus beschrieben in den obigen Abschnitten (wobei nur EINE
Aufnahme gemacht wird) erlaubt es nicht dass alle 3D Kalibrierpunkte
in der gleichen Fläche liegen.
Zhang 2000, und auch die Beschreibung in Wöhler 2009 (p. 22) eines
Papers von Tsai aus 1987, benutzen jedoch planare calibration rigs.
Sie nehmen das Welt-Koordinatensystem W so dass das Kalibriermuster in
der Fläche zW = 0 liegt. D.h. by construction hat man zW = 0 für alle
3D Kalibrierungspunkte.
Gleichung M4 wird
+- -+ +- -+
| xS~ | | xW |
| yS~ | = M | yW |
| zS~ | | 0 |
+- -+ | 1 |
+- -+
Die dritte Spalte von M "fällt weg".
Gleichung M5 (weiter oben), mit den ausgeschriebenen Matrixelementen,
reduziert sich zu
+- -+ +- -+ +- -+ +- -+
| xS~ | | -fx 0 Ox | | r11 r12 -tRx | | xW |
| yS~ | = | 0 -fy Oy | | r21 r22 -tRy | | yW |
| rS~ | | 0 0 1 | | r31 r32 -tRz | | 1 |
+- -+ +- -+ +- -+ +- -+
\____________/ \_______________/
Mint Mext
\___________________________________/
M
(Gleichung Z1)
wo die dritte Spalte (r13,r23,r33) in Mext weggefallen ist. Matrix M =
Mint Mext reduziert sich zu einer 3x3 Matrix.
Da die Gleichungen M1a, M1b hier noch immer gelten, ist auch die 3x3
Matrix M nur definiert up to an arbitrary scaling factor. Ähnlich wie
oben für den nicht-planaren Fall kann man aus einer bekannten Menge
Korrespondenzen (aus einer Aufnahme des calibration rigs) die 9 Elemente
m_ij dieser 3x3 Matrix zu feststellen up to a common multiplicative
factor (Vorgang ist wieder über A m = 0 und SVD, mit danach eventuell
einem Levenberg-Marquard refinement, alles wie oben beschrieben).
Diese 9 Zahlen, mit 8 Freiheitsgraden, reichen auf jedem Fall aus
dazu, aus der single Aufnahme ein Kameramodell zu kalibrieren das 8
parameters hat.
Mehrere der gelesenen Literaturstücke schreiben dass Ox = Oy = 0
für viele Kameras eine gute erste Annäherung ist. Das verringert die
Anzahl intrinsic parameters zu 2, nämlich { fx, fy }, also mit den
extrinsic parameters insgesamt 8 Kamera-parameters. Dieses Kameramodell
würde man kalibrieren können aus einer einzigen Aufname eines planaren
Kalibriermusters.
[ Bemerkung: Ich vermute, dass ein so kalibriertes Modell nicht
genau sein würde für 3D scene points die sich befinden auf einer zW weit
entfernt von 0, d.h. weit entfernt von der Fläche des planar calibration
rigs. ]
Zhang 2000 und Wöhler 2009 (p. 22) beschreiben Verfahren wie man mit
einem planaren calibration rig trotzdem ein Kameramodell kalibrieren
kann das MEHR als 2 intrinsic parameters hat. Alle dieser Verfahren
beruhen (nicht erstaunlicherweise) darauf, dass man MEHRERE Aufnahmen
macht, und den planaren calibration rig zwischen Aufnahmen VERSCHIEBT
und/oder DREHT relativ zur Kamera.
Sowohl Zhang 2000 wie auch Wöhlen 2009 (p. 22) nutzen aus dass
wir in Gleichung Z1 ein 'a priori' Wissen haben dass in den obigen
Abschnitten noch nicht benutzt wurde; nämlich wir wissen dass die Spalten
(r11,r21,r31) und (r12,r22,r32) in Mext, Spalten aus der Rotationsmatrix
R, orthonormal sind.
Unten beschreibe ich den Algorithmus von Zhang 2000. Die Beschreibung
erklärt auch wieviele Aufnahmen benötigt sind, und warum der planare
calibration rig zwischen Aufnahmen verschoben/rotiert werden muss.
>>> Zhang -- Übersicht
Zhang 2000 benutzt eine planare calibration rig, wovon er n Aufnahmen
macht, unter unterschiedlichen Rotationen u/o Translationen relativ zur
Kamera. Die Aufnahmen werden alle mit der selben Kamera gemacht,
nämlich die Kamera wovon die intrinsic parameters festzustellen sind.
Meine Eindruck ist, dass er sein Verfahren vor allem abgezielt hat auf
das feststellen der _intrinsic_ parameters einer Kamera. (Sein Verfahren
kann bis zu 5 intrinsic parameters feststellen.)
Jede Aufnahme betrachtet er zunächst separat. Folgende sind die Schritte
die er durchführt für jede Aufname separat:
Er nimmt das W-Koordinatensystem so dass zW = 0 die Fläche seines
planaren calibration rigs ist, und benutzt die "reduzierte"
Vektoren und Matrizen aus Gleichung Z1. Er findet eine Menge
Korrespondenzen zwischen (xW,yW) auf seinem calibration rig und
Bild-Punkten (xS,yS). Aus diese Korrespondenzen stellt er die
Werte fest der Elemente m_ij in der 3x3 Matrix M (in Gleichung
Z1), up to a common multiplicative factor.
Bemerke dass jede Aufname ein anderes W-Koordinatensystem benutzt,
da der calibration rig zwischen aufnahmen verschoben u/o gedreht wird.
Zhang benutzt keine Information darüber wie die W-Koordinatensysteme
der unterschiedlichen Aufnahmen sich zu einander beziehen.
Die m_ij die so gefunden werden für jede separate Aufnahme sind eine
Schätzung (Messung) für die Matrix M = Mint Mext die gültig ist für
die _jeweilige_ Aufnahme. Notation: *M = matrix der geschätzten
(gemessenen) m_ij für eine Aufnahme.
Menge Korrespondenzen --> *M für Aufnahme 1
aus Aufnahme 1
Menge Korrespondenzen --> *M für Aufnahme 2
aus Aufnahme 2
... ...
Zwischen Aufnahmen ändert sich Mext. Zhang nimmt an, dass sich zwischen
Aufnahmen Mint nicht ändert. (Das bedeutet m.E. dass die Aufnahmen nicht
lange Zeit nach einander genommen werden müssen, so dass Faktoren die
die intrinsic parameters beeinflüssen könnten, wie z.B. die Temperatur,
sich zwischen Aufnahmen nicht stark ändern.)
Die Situation ist also die folgende:
(Subscript _i angehängt an Größe die gilt für Aufnahme i)
Menge Korrespondenzen --> *M_1 = Mint Mext_1
aus Aufnahme 1
Menge Korrespondenzen --> *M_2 = Mint Mext_2
aus Aufnahme 2
... ...
wo Mint für jede Aufnahme die gleiche ist. Alle *M_i sind bekannt.
Mint, und alle Mext_i, sind unbekannt. Bemerke weiter dass jede der
3x3 matrices Mext_i ein anderes W-Koordinatensystem benutzt.
Die Idee hinter Zhang's Methode ist nun, diese n Matrixgleichungen
*M_1 = Mint Mext_1
...
*M_n = Mint Mext_n
zu behandeln als n "Messpunkte" für Mint. Damit dies möglich ist,
sollte aus jedem dieser n "Messpunkte" die matrix Mext_i irgendwie
eliminiert werden. Das gelingt dadurch, dass wir ausnutzen dass wir
wissen dass in jeder Mext_i die ersten zwei Spalten (r11,r21,r31) und
(r12,r22,r32) orthonormal sind.
>>> Zhang -- Elimination von Mext_i aus *M_i = Mint Mext_i
Innerhalb dieser subsection betrachte ich nur _eine_ Aufnahme separat.
Dass die m_ij bekannt sind up to a common multiplicative factor heißt:
lambda *M = Mint Mext
(wo *M = die 3x3-matrix der m_ij für die betrachtete Aufnahme, und
Mext = die 3x3-matrix Mext für die betrachtete Aufnahme.)
Ausschreiben in Spalten ergibt:
lambda [ m1 m2 m3 ] = Mint [ r1 r2 tR ]
(Gleichung Z2)
wo
mi = i-te Spalte aus matrix *M
ri = i-te Spalte aus der rotation matrix R für die betrachtete Aufnahme
tR = translation vector für die betrachtete Aufnahme.
Das Wissen das wir jetzt ausnutzen möchten ist dass die Vektoren r1
und r2 orthonormal sind:
r1^T r2 = 0 (Gleichung Z3a)
r1^T r1 = r2^T r2 = 1 (Gleichung Z3b)
(Notation: T = transpose, ^ = superscript.)
Aus Gleichung Z2 entnehmen der Spalten für r1 und r2 ergibt:
lambda m1 = Mint r1
lambda m2 = Mint r2
Die Spalte tR aus Mext, und die Spalte m3 aus *M, werden für das
Feststellen von Mint nicht benutzt. Die matrix Mint (mit den zu
bestimenden intrinsic parameters) ist klein, und analytsisch umrechnen
zwischen Mint und seine Inverse Mint^-1 ist einfach. Wir dürfen daher
úmschreiben nach
r1 = lambda Mint^-1 m1 (Gleichung Z4a)
r2 = lambda Mint^-1 m2 (Gleichung Z4b)
Hiermit können wir das Scalarprodukt von ri und rj schreiben als
ri^T rj
= lambda^2 ( Mint^-1 mi )^T ( Mint^-1 mj )
= lambda^2 mi^T Mint^-T Mint^-1 mj
Substitution unserer r1 und r2 aus Gleichungen Z4a, Z4b in die
Orthogonalitätsbedingungen (Gleichungen Z3a, Z3b) ergibt also:
m1^T Mint^-T Mint^-1 m2 = 0 (Gleichung Z5a)
m1^T Mint^-T Mint^-1 m1 = m2^T Mint^-T Mint^-1 m2
(Gleichung Z5b)
Das Ziel ist erreicht. Z5a, Z5b sind Gleichungen in denen nur noch
vorkommen: unsere unbekannte Mint und die bekannten Messwerte m1, m2
(= Spalten aus *M) für die betrachtete Aufnahme.
In unseren Gleichungen Z5a, Z5b kommt die Kombination ( Mint^-T Mint^-1 )
ständig wiederholt vor. Dieses Matrixprodukt ist eine symmetrische 3x3
matrix die nur von den intrinsic parameters abhängt. [ Symmetrisch weil
(A^T A)^T = A^T A für jede Matrix A. ] Auch das Umrechnen zwischen
Mint^-1 und der Produkmatrix ( Mint^-T Mint^-1 ) ist analytisch einfach
machbar. Wir führen daher ein die Notation
B = Mint^-T Mint^-1
so dass wir Z5a, Z5b schreiben können als
m1^T B m2 = 0 (Gleichung Z6a)
m1^T B m1 - m2^T B m2 = 0 (Gleichung Z6b)
Da die 3x3 matrix B symmetrisch ist:
+- -+
| b11 b12 b13 |
B = | b12 b22 b23 |
| b13 b23 b33 |
+- -+
enthalten die Gleichungen Z6a, Z6b nur 6 unknowns: b11, b12, b13, b22,
b23, b33.
>>> Zhang -- Ableiten von Mint aus den *M_i für die n Aufnahmen
Jede Aufnahme gibt uns ein Gleichungspaar Z6a, Z6b. Es werden n
Aufnahmen gemacht, so dass wir n solcher Gleichungspaare erhalten:
m1_1^T B m2_1 = 0
m1_2^T B m1_1 - m2_1^T B m2_1 = 0
m1_2^T B m2_2 = 0
m1_2^T B m1_2 - m2_2^T B m2_2 = 0
...
m1_n^T B m2_n = 0
m1_n^T B m1_n - m2_n^T B m2_n = 0
(Gleichungssystem Z7)
wo
m1_i = erste Spalte aus *M für Abbildung i
m2_i = zweite Spalte aus *M für Abbildung i.
Die matrix B (die nur von den intrinsic parameters abhängt) ist gleich
für jede Abbildung. Unsere unknowns sind die sechs Elemente b11, b12,
b13, b22, b23, b33 der Matrix B. Also insgesamt 2n Gleichungen und
6 unknowns.
Ein einfaches Ausschreiben dieser Gleichungen in den Elementen von B,
m1_i und m2_i zeigt dass die Gleichungen in obiges Gleichungssystem sind
alle linear in den unknowns b11, b12, b13, b22, b23, b33. (Und dass
die Koëfficienten der b_ij Summen sind von Produkten der Elemente von
m1_i und m2_i.)
Gleichungssystem Z7 kann daher geschrieben werden als Matrixgleichung
C b = 0 (Gleichung Z8)
wo
b = Vektor (b11, b12, b13, b22, b23, b33) mit unseren unknowns
C = 2n x 6 matrix mit bekannten Koëfficienten.
Das genau wie oben im Abschnitt
"Kalibrierung aus N Korrespondenzen, ohne radial distortion"
wieder eine homogene Matrixgleichung. Genau wie dort finden wir wieder
den optimalen b aus singular value decomposition C = U D V^T. Der b
der den error || e || = || C b || minimiert ist der Spaltenvektor v1
aus V der gehört zur kleinsten singular value s1 (aus D).
Die matrix C hat (theoretisch) höchstens den rank 5 (nullspace
1-dimensional). Damit C diesen maximalen rank erreicht wird, braucht
man mindestens 5 Gleichungen, das heißt (mit 2 Gleichungen pro Aufnahme)
mindestens 3 Aufnahmen. Die 6 Zahlen b_ij sind damit festgelegt up to
an arbitrary common factor. Das reicht aus dazu, 5 intrinsic parameters
festzulegen (zu kalibrieren).
[Siehe für Kameramodelle mit weniger als 5 intrinsic parameters die
subsection "Zhang -- Spezialfall für 4 intrinsic parameters" unten.]
>>> Zhang -- Verschieben/drehen des calibration rigs zwischen Aufnahmen
Wie immer bei homogenen linearen Gleichungssystemen ist es wichtig dass
die nullspace von C in Gleichung Z8 nicht mehr Dimensionen haben sollte
als 1, d.h. nur ein singular value hat gleich 0. [Wegen Messfehler
erwarten wir auch hier wieder dass es ein singular value gibt nahe 0,
mit allen anderen singular values deutlich > 0.]
Die Berechnungen sollten m.E. auch hier wieder die singular values
überprüfen, und Alarm geben wenn die singular values nicht sind wie
benötigt (einer near 0, alle andere deutlich > 0).
In der Matrix C werden mehr als ein singular value (nahe) null-wertig
wenn die m1_i und m2_i für unterschiedliche Aufnahmen (nahezu) gleich
sind.
m1_i und m2_i sind Spalten aus der "gemessenen" matrix *M_i für
die Aufnahme i, die abgeleitet wurde aus der Menge Korrespondenzen
gefunden in Aufnahme i. Wir wissen dass
*M_i = Mint Mext_i
wo Mint konstant ist für alle Aufnahmen. Die *M_i der unterschiedlichen
Aufnahmen sind daher nur dann unterschiedlich, wenn die Mext_i --
die rotation und/oder die translation des calibration rigs relativ zur
Kamera -- in den Aufnahmen unterschiedlich sind.
Beispiel: Wenn in allen n Aufnahmen der calibration rig exact die
gleichen Position und Orientierung relativ zur Kamera hat -- oder
wenn man n mal die gleiche Aufnahme benutzt --, dann ist Mext_i
für jede Aufnahme gleich, und wir erhalten für jede Aufnahme die
gleiche *M_i. Gleichungssystem Z7 besteht jetzt aus eine n-fache
Wiederholung von zwei Gleichungen; und Matrix C in Gleichung z8 hat
jetzt denselben rank wie eine C die nur die zwei Gleichungen für
eine einzige Aufnahme enthält, nämlich rank = 1.
Obiges Beispiel ist ein Extremfall. Wenn man das Kalibriermuster
zwischen den Aufnahmen nur ein sehr wenig verschiebt, dann wird die
matrix C immer noch "nahezu" den rank 1 haben (d.h. nur eine singular
value haben der deutlich > 0 ist).
Die Mindest-Bedingung die erfüllt werden muss für die Größe der
Verschiebung/Rotation des Kalibriermusters zwischen Aufnahmen ist
m.E. dass C nur ein singular value nahe 0 hat, und dass alle andere
singular values eine Größe haben die "über die Mess-Ungenauigkeit liegt"
-- die "signal to noise ratio" sollte hoch genug sein. Offene Fragen:
Wie genau definiert und schätzt man die Mess-Ungenauigkeiten ?, und
wie bringt man man die Mess-Ungenauigkeiten und die singular values auf
der gleichen scale damit man sie vergleichen kann ? Fürs Lösen dieser
Fragen würde ich mehr Zeit brauchen.
Die Qualität (Genauigkeit) der erhaltenen Kamera-parameters verbessert
sich wenn die "signal to noise ratio" sich steigert. Ich denke (vermute
sehr stark) dass für jede Kalibration gilt dass die 3D Kalibrierpunkte
am besten gleichmäßig verteilt sein sollten über den ganzen 3D Raum
worin die 3D scene of interest sich befindet; so würde man die maximale
Genauigkeit in den Kamera-parameters erreichen.
Daraus schließe ich dass die Verschiebungen u/o Drehungen des planaren
Kalibriermusters zwischen den Aufnahmen am besten so sein sollten
dass die 3D Kalibrierpunkte, insgesamt über allen Aufnahmen, sich weit
verteilen über den 3D-Raum of interest.
>>> Zhang -- Spezialfall für 4 intrinsic parameters
Zhang 2000 zeigt dass man die 4 intrinsic parameters aus dem Kameramodell
+- -+
| -fx 0 Ox |
Mint = | 0 -fy Oy |
| 0 0 1 |
+- -+
aus nur ZWEI Aufnahmen kalibrieren kann. Analytisch ableiten der
Matrix B = ( Mint^-T Mint^-1 ) aus dem obigen mInt ergibt Ausdrücke
für die Elemente von B in terms of den 4 intrinsic parameters. Mit der
obigen Mint ergibt das b12 = 0, unabhängig von den intrinsic parameters.
Das heißt dass "theoretisch" der unknown b12 null sein sollte.
Mit anderen Worten, wir können die Gleichung
b12 = 0
hinzufügen in unser Gleichungssystem Z7. Das so erweiterte
Gleichungssystem Z7 hat dann bei 2 Aufnahmen insgesamt 5 Gleichungen.
Matrix C wird eine 5 x 6 matrix. Da C b = 0, hat C theoretich maximal
den rank 4. (Dieser rank wird nur erreicht wenn in den 2 Aufnahmen die
calibration rig unterschiedlich gedreht/positioniert wird.) Das reicht
dazu aus, die 4 intrinsic parameters zu bestimmen.
[Bemerkung: EINE einzige Aufnahme eines planaren Kalibriermusters
reicht aus für 2 intrinsic parameters, aber dafür braucht man das
Zhang-Verfahren nicht -- siehe den Anfang dieses Abschnittes "Planarer
calibration rig".]
>>> Zhang -- Extrinsic parameters
Da der calibration rig zwischen Aufnahmen verschoben/gedreht wird,
sind die extrinsic parameters für jede Aufnahme unterschiedlich.
Die extrinsic parameters für Aufnahme i sind die rotation R und
translation tR des W-Koordinatensystems -- das im calibration rig liegt
-- relativ zum Kamera-Koordinatensystem.
Gegeben die intrinsic parameters, kann man diese R und tR für jede
Aufnahme berechnen wie folgt. Hierzu wird die Information über r1,
r2, tR aus Gleichung Z2 benutzt die bei der Bestimmung der intrinsic
parameters ignoriert wurde.
r1/lambda und r2/lambda erhält man aus Gleichungen Z4a und Z4b.
Normieren von r1 und r2 nach unit length ergibt den scale factor lambda.
r3 = r1 x r2. Aus Gleichung Z2 entnehmen der tR-Spalte ergibt:
lambda m3 = Mint tR
woraus tR bekannt ist.
>>> Zhang -- Verfeinerung
Genauso wie oben beschrieben, im Abschnitt "Verfeinerung der
Kamera-parameters" für den nicht-planaren Fall, behandelt auch Zhang
die auf diese weise "analytisch" erhaltenen Kamera-parameters nur als
initial guess für ein Levenberg-Marquard refinement.
Der Grund dafür ist auch bei Zhang genau der gleiche wie vorher, nämlich:
(1) der Fehler || e || = || C b || der minimiert wird in der analytischen
Methode ist eine arbiträre, nicht physisch sinvolle Größe;
(2) die analytische Methode erlaubt keine Schätzung der radial distortion
parameters.
Das von Zhang benutzte nicht-linearen Levenberg-Marquard refinement
geht genau wie im von mir oben beschriebenen nicht-planaren Fall:
- Er benutzt während des refinements das volle (nicht linearisierte)
Kameramodell mit radial distortion (initial guess wieder: k1 = k2 = 0);
- Während des refinements ist der error der minimiert wird wieder
E = SUM_i [ d_i^2 ] genau wie in Gleichung K3 (im Abschnitt
"Verfeinerung der Kamera-parameters" oben).
---------------------------------------------------------------------------
Anmerkungen
-----------
[note 1] Homogene Koordinaten:
Hier kommt mir die Idee, über einen Trick ähnlich wie homogene
Koordinaten auch die radial distortion Formeln zu "linearisieren".
Vielleicht so:
+- -+ +- -+ +- -+
| xPV | = | alfa xPN | <----> | xPN |
| xPV | | alfa yPN | | yPN |
+- -+ +- -+ | alfa |
+- -+
+- -+ +- -+ +- -+
| xPN | | 1 0 0 0 0 | | xPN |
| yPN | = | 0 1 0 0 0 | | yPN |
| alfa | | 0 0 1 k1 k2 | | 1 |
+- -+ +- -+ | r^2 |
| r^4 |
+- -+
In Wöhler, Abschnitt 1.4.2. "The Direct Linear Transform Method"
(p. 20) scheint eine Form dieser Idee angewendet zu werden.
Ich habe leider noch keine Zeit gefunden, den Nutzen und
Anwendbarkeit dieser Idee durchzuanalysieren.
[note 2] Braucht 6 Korrespondenzen:
Dieser Kalibrier-Algorithmus braucht also eine Korrespondenz mehr
als die 5 Korrespondenzen die am Anfang des Kapitels "KALIBRIERUNG"
vorhergesagt wurden als das Minimum fürs Festlegen unserer 10
Kamera-parameters. Die 11 Freiheisgrade für die m_ij die wir in
dem von mir hier beschriebenen Algorithmus alle füllen, werden
reduziert zu den 10 unserer Kamera-parameters während wir am Ende
des Abschnittes "Kalibrierung aus N Korrespondenzen, ohne radial
distortion" die Kamera-parameters ableiten aus den m_ij und dabei
die Rotationsmatrix R normalisieren.