Camera🚧

Camera

3D物體轉換成2D圖片的機制。

甚至從2D圖片推敲3D物體。

Camera的領域學科

數位相機的不同部件源自不同領域。

鏡頭:光學工程之optical design
成像:光電工程之image sensor
攝影:交叉學科computational photography
測量:土木工程/空間資訊工程之photogrammetry、計算機科學之3D vision

3D Vision

本文只針對3D vision。流程大致如下:

我覺得應該先來一張精美的流程圖。有空再來畫。

1. image feature extraction
   一張相片,找到引人注目的地方,求得確切位置。

2. image feature description
   一張相片,引人注目的地方,重新表示成數值,方便比對。

3. image feature matching
   兩張相片,引人注目的地方,求得對應關係。

4. camera calibration
   事先測量物體表面座標。
   一臺相機,拍一張相片,求得相機的位置、方向、焦距。

5. camera motion estimation
   一臺相機一直動,視作多臺相機。
   多臺相機,各拍一張相片,求得每臺相機的位置、方向、焦距。

6. surface point cloud generation
   兩張相片,找到物體表面座標。

7. surface point cloud registration
   多張相片,對齊物體表面座標。

8. surface reconstruction
   物體表面座標重新拼成網格/參數曲面/貝茲曲面/體素/等值曲面/距離場。
1. structure from motion
   一臺相機一直動,得到場景模型。即是上述八項步驟。

2. match move
   一臺相機一直動,把東西塞進去場景裡面,配合相機方位。

3. augmented reality
   塞進去的東西可以跟場景互動,例如碰撞偵測。

使用現成軟體進行攝影測量

資料集:KITTIMiddleburyETH3D

學術專案:ORB-SLAMOpenSfMColmap

函式庫:OpenCVBoofCVjsfeatscikit-image

工具:MATLAB Computer Vision ToolboxKorniaARCoreUnity

軟體:Blender3DEqualizerSyntheyes

公司:相機KODAK、成像設備Canon、手機Apple。

產業:相機模組製造商航空攝影公司視覺特效公司

Image Tracking🚧

簡介

Image Keypoint Detection

要點偵測。擷取圖片引人注目的地方。

此處的圖片是臺北郵局。我為了排版美觀,事先縮小圖片尺寸。

Hessian keypoint detector:每個像素求得det(Hessian matrix),擷取數值偏高者。這些像素很可能是陡變突顯之處。

Hessian matrix是二次微分。Hessian matrix的兩個特徵值大致是主曲率。Hessian matrix的兩個特徵值的乘積大致是高斯曲率。因為函數曲面和參數曲面的曲率公式不一樣,所以說大致。

主曲率可以判斷凹凸。兩個特徵值λ₁和λ₂,皆正隆起(有局部極大值)、皆負陷落(有局部極小值)、一正一負馬鞍面、一零折彎面、皆零斜平面。

高斯曲率可以判斷凹凸。兩個特徵值的乘積det(Hessian matrix) = λ₁λ₂,數值偏高者,隆起陷落明顯。

Laplacian keypoint detector:改為tr(Hessian matrix),擷取數值偏高者與偏低者。亦可取絕對值或取平方,擷取數值偏高者。

tr(Hessian matrix) = λ₁+λ₂是Laplacian。Laplacian是梯度的散度。Laplacian大致是主曲率的總和、平均曲率的兩倍。

平均曲率可以判斷曲直。數值偏高者,彎曲明顯。

換句話說:套用Laplacian filter,擷取數值偏高者與偏低者。

     ∂I
Iₓ = ——                          first-order derivative
     ∂x

       ∂²I
Iₓ₝ = —————                      second-order derivative
      ∂x ∂y

I(x,y)                           R/G/B value at pixel (x,y)

Iₓ(x,y) = I(x,y) - I(x-1,y)      x gradient
I₝(x,y) = I(x,y) - I(x,y-1)      y gradient

Iₓ₝(x,y) = Iₓ(x,y) - Iₓ(x,y-1)   y gradient of x gradient
I₝ₓ(x,y) = I₝(x,y) - I₝(x-1,y)   x gradient of y gradient

H(x,y) = ⎡ Iₓₓ(x,y) I₝ₓ(x,y) ⎤   Hessian matrix
         ⎣ Iₓ₝(x,y) I₝₝(x,y) ⎦   (second-order derivative matrix)

H = ⎡ Iₓₓ I₝ₓ ⎤                  skip (x,y)
    ⎣ Iₓ₝ I₝₝ ⎦

det(H) = Iₓₓ I₝₝ - Iₓ₝ I₝ₓ       determinant of Hessian matrix
       = Iₓₓ I₝₝ - (Iₓ₝)²        since Iₓ₝ = I₝ₓ

tr(H) = Iₓₓ + I₝₝                trace of Hessian matrix
                                 (Laplacian)

Image Blob Detection

斑點偵測。引人注目的地方,從像素拓寬成區域。

LoG blob detector:斑點形狀是自訂圓形。Laplacian keypoint detector拓寬成圓形。一、指定圓形半徑。二、針對一種圓形中心位置,擷取圓形裡面的像素。三、每個像素求得Laplacian。四、套用窗函數Gaussian window,中央高外圍低,專注於圓形中心附近的Laplacian。五、求得Laplacian總和。六、窮舉圓形中心位置。當總和是局部極大值、局部極小值,推定為斑點。

換句話說:套用Gaussian of Laplacian filter(等同Laplacian of Gaussian filter,卷積交換律)。自訂濾波器寬度(變異數),推定為斑點半徑。新圖片局部極值所在之處,推定為斑點中心位置。

援引高斯散度定理:圓形裡面Laplacian總和=圓形邊界梯度投影到法向量之總和。此演算法本質是梯度通量偵測,而且圖片事先平滑化。也有人解釋為台地盆地偵測。

DoG blob detector:LoG blob detector的高速近似算法。

Kadir–Brady saliency detector:斑點形狀是自訂造型。窮舉斑點位置。計算斑點的histogram與entropy,擷取entropy較大者。

maximally stable extremal regions (MSER):自動尋找斑點形狀。枚舉臨界值,分別二值化。追蹤每個connected component的面積,擷取面積變化緩慢者。

Image Blob Adaptation

斑點適應。引人注目的地方,調整尺寸和形狀。

automatic scale selection:調整blob尺寸。一、blob改成各種尺寸。二、計算blob裡面所有像素的梯度長度的平均數。三、平均數依照blob尺寸排列,形成一串數列。四、局部極大值所在之處,表示梯度變化最大,推定為正確尺寸。當局部極大值不只一個,推定blob不存在。

affine shape adaptation:調整blob形狀,成為橢圓形。取得blob裡面所有像素,求得二次動差矩陣總和,實施特徵分解,得到橢圓形,作為新形狀。重複上述步驟,直至blob形狀不變。

換句話說:一個像素的梯度,視作一個二維點座標。blob裡面所有像素的梯度,視作一群二維點座標。principal component analysis(不必預先減去平均數),得到橢圓形。

I(x,y)                           R/G/B value at pixel (x,y)

Iₓ(x,y) = I(x,y) - I(x-1,y)      x gradient
I₝(x,y) = I(x,y) - I(x,y-1)      y gradient

P = ⎡ Iₓ(x₁,y₁) Iₓ(x₂,y₂) ... Iₓ(xₙ,yₙ) ⎤
    ⎣ I₝(x₁,y₁) I₝(x₂,y₂) ... I₝(xₙ,yₙ) ⎦

PPᵀ = ⎡ ∑IₓIₓ ∑I₝Iₓ ⎤ = ∑ S(x,y)
      ⎣ ∑IₓI₝ ∑I₝I₝ ⎦

S(x,y) = ⎡ IₓIₓ I₝Iₓ ⎤           second moment matrix
         ⎣ IₓI₝ I₝I₝ ⎦           (structure tensor)

詳細原理請見這篇講義。一個像素與另一個鄰近像素的數值差異,採用平方誤差,可以改寫成二次型。推廣成區域,仍然可以改寫成二次型。二次型的最大值/最小值,位於最大/最小的特徵值的特徵向量。以特徵向量作為橢圓形的主軸/次軸,以特徵值比例作為橢圓形的軸長比例,藉此得到差異最大的區域、引人注目的地方。

  E(x,y)                                       error

= (I(x+Δx,y+Δy) - I(x,y))²                     squared error

≈ (I(x,y) + Iₓ(x,y)Δx + I₝(x,y)Δy - I(x,y))²   first-order
                                               Taylor approximation
= (I + IₓΔx + I₝Δy - I)²                       skip (x,y)

= (IₓΔx + I₝Δy)²                               elimination

= (IₓΔx)² + (I₝Δy)² + 2(IₓΔx)(I₝Δy)            expansion

= [ Δx Δy ] ⎡ IₓIₓ IₓI₝ ⎤ ⎡ Δx ⎤               quadratic form
            ⎣ IₓI₝ I₝I₝ ⎦ ⎣ Δy ⎦

= [ Δx Δy ] S(x,y) ⎡ Δx ⎤                      second moment
                   ⎣ Δy ⎦                      matrix

註:Δx是物理學displacement,不是數學finite difference!
    ∑   E(x,y)                                 sum of errors
 (x,y)∈R                                       in the region

≈   ∑    [ Δx Δy ] ⎡ IₓIₓ IₓI₝ ⎤ ⎡ Δx ⎤        substitution
 (x,y)∈R           ⎣ IₓI₝ I₝I₝ ⎦ ⎣ Δy ⎦

= [ Δx Δy ] ⎛    ∑    ⎡ IₓIₓ IₓI₝ ⎤ ⎞ ⎡ Δx ⎤   linearity
            ⎝ (x,y)∈R ⎣ IₓI₝ I₝I₝ ⎦ ⎠ ⎣ Δy ⎦

= [ Δx Δy ] ⎛    ∑   S(x,y) ⎞ ⎡ Δx ⎤           sum of second
            ⎝ (x,y)∈R       ⎠ ⎣ Δy ⎦           moment matrix
  argmin    ∑   E(x,y)   subject to (Δx)² + (Δy)² = 1
   Δx,Δy (x,y)∈R

= eigvec    ∑   S(x,y)
         (x,y)∈R

圖片的數學知識:X-Invariant Transformation

我們希望圖片經過位移、縮放、旋轉,keypoint和blob仍然相同。分別稱作位移不變、縮放不變、旋轉不變。

設計演算法的時候,採用特殊手法,可以達成位移不變、縮放不變、旋轉不變。以下介紹這些手法。

translation-invariant:只使用鄰近像素。鄰居關係不因位移而改變。凡是image filter都是位移不變。

scale-invariant:圖片縮放為各種尺寸,從中找到最佳尺寸。如此一來,計算結果不因縮放而改變。

rotation-invariant:圖片旋轉為各種角度,從中找到最佳角度。如此一來,計算結果不因旋轉而改變。

圖片的數學知識:X-Invariant Transformation有兩種定義

X-invariant transformation,嚴格來說擁有兩種定義,視情況而定。為了行文流暢,大家一般不會挑明是哪種定義。

transformation: F(x) = y
definition 1:   F(X(x)) = y      輸入套用X,輸出仍然相同。
definition 2:   F(X(x)) = X(y)   輸入套用X,等同於輸出套用X。

先前小節提到的演算法,有些已經達成位移不變、縮放不變、旋轉不變。以下進行詳細介紹。

keypoint:針對連續函數,det(Hessian matrix)tr(Hessian matrix)恰是旋轉不變(定義一)。針對離散陣列,事情相當複雜,必須額外考慮image warping的演算法。於是大家轉頭看向連續函數,姑且當作是旋轉不變(定義一)。

blob:LoG blob detector是位移不變暨旋轉不變(定義二)。automatic scale selection追加了縮放不變(定義二)。affine shape adaptation將圓形拓寬為橢圓形,同時維持三種不變。

affine shape adaptation:二次動差矩陣總和,實施特徵分解,特徵值是旋轉不變(定義一),特徵向量是旋轉不變(定義二)。

Image Similarity-Invariant Blob Detection

相似不變斑點偵測=斑點偵測+斑點適應。

位移不變、縮放不變、旋轉不變,同時達成,合稱相似不變。

Harris affine region detector:達成相似不變。採用下個小節的Harris corner detector。圖片縮放成各種尺寸,每張圖都找corner。每一個corner,以自己為中心,嘗試找橢圓形blob。圓形拓寬為橢圓形,同時維持三種不變。

此處的affine是仿射變換(圓形拓寬為橢圓形),不是仿射不變(圖片實施仿射變換,橢圓形blob仍然相同)。

maximally stable extremal regions (MSER):涵蓋相似不變。connected component具備各種不變,包括位移縮放旋轉歪斜投影鏡射,包括剛體相似線性仿射射影拓樸。包山包海反而抓不到重點。

center surround extremas (CenSurE):放棄旋轉不變。放棄圖片縮放,採用區域縮放。放棄濾波器,採用區域總和。放棄梯度,採用原始數值。放棄圓形,採用特殊區域造型,例如兩個正方形構成的八角星。一切都是為了算得更快。

圖片的數學知識:大家需要的X-Invariant Transformation

實際需求如下。目前沒有經典演算法。太麻煩了,不如不做。

projective-invariant:射影不變。相機位置導致交角變化。

distortion-invariant:扭曲不變。鏡頭造型導致彎曲變化。

illumination-invariant:照明不變。光源位置導致色彩變化。

Image Corner Detection

角點偵測。找到圖片景物的轉角。

Harris corner detectorShi–Tomasi corner detector:利用affine shape adaptation的特徵值大小來判斷轉角。一、指定區域造型。二、針對一種區域中心位置,擷取區域裡面的像素。三、每個像素求得梯度。四、套用窗函數Gaussian window,專注於區域中心附近的梯度。五、求得二次動差矩陣總和,實施特徵分解。即是梯度們的principal component analysis(不必預先減去平均數)。六、以兩個特徵值的絕對值大小,判斷該區域是否有角。兩值都足夠大,就表示梯度們有兩個主要方向,就表示有兩條較直的邊緣,就表示有角。七、判斷條件可以改寫成determinant和trace。

features from accelerated segment test (FAST):一、待測像素,作為圓心。像素值加上/減去特定數值,作為臨界值。二、擷取圓周16點,若連續12點大於/小於臨界值,視作角。三、預先擷取圓周上下左右4點,若連續3點大於/小於臨界值,才做二,節省時間。

Image Edge Detection

邊緣偵測。找到圖片景物的邊緣。

Marr–Hildreth edge detector:圖片套用LoG filter或者DoG filter,過零處(正負變化之處)當作邊緣。

Canny edge detector:一、平滑兼去噪:Gaussian filter。二、求邊緣:求得梯度長度、梯度角度。三、找到最明顯的邊緣:根據梯度角度,取得前後兩個像素。比較梯度長度,若是局部極大值(前後皆較小),則保留,否則設為零。四、清除不明顯的邊緣:設定門檻,梯度長度太短者設為零。五、銜接邊緣。

Image Line Segment Detection

線段偵測。找到圖片景物的平整邊緣。

Hough line transform:偵測直線。一、求邊緣:求得梯度長度,設定門檻。二、針對一個像素(x,y),窮舉穿過此像素的直線們,求得原點在直線上的投影點,再換成極座標表示法(r,θ),一條直線對應一個點,直線們對應一條曲線r = xcosθ + ysinθ。一個像素得到一條曲線。三、圖片每個像素皆換成曲線。曲線們的交點,即是穿過大量像素的直線!四、建立細膩方格紙(離散化),曲線穿過方格,則投票加一,得票多的方格當作交點。五、轉換過程是三步驟:對偶、轉成極座標、投票。

Hough cycle transform:偵測圓。窮舉各種半徑。針對一個像素(x,y),窮舉穿過此像素的所有圓;圓穿過像素,則投票加一,得票多的像素就是圓心。

line segment detector (LSD):偵測直線線段。找到梯度長度大於臨界值的像素。連接梯度角度差異小於臨界值的像素。形成矩形外框,框住直線線段。當像素間距大於臨界值,則拆散像素,重來一遍,形成更小矩形外框。最後以原創公式NFA篩選矩形。

edge drawing (EDLines):LSD改良版本。放棄矩形外框,採用直線擬合。誤差大於一的像素,當作不在直線上。

Image Feature Detection

特徵偵測。找到圖片的鋩鋩角角。

keypoint      要點/關鍵點
blob          斑點
corner        轉角/角點
edge          邊緣
ridge         稜線/脊
line segment  線段

Image Feature Extraction

特徵萃取。挑選容易認知的鋩鋩角角。

大家將此主題併入image feature description。請見下個小節。

也有人認為此主題等同image feature engineering。也就是說,本章所有小節,依序跑過一遍。

histograms of oriented gradients (HOG)的作法:我全都要。圖片分解成大量16×16正方形區域,交疊一半。

scale-invariant feature transform (SIFT)的作法:Harris affine region detector得到大量橢圓形blob。一個橢圓形實施旋轉縮放,擺正橢圓形,形成16×16正方形區域。

然而橢圓形blob不適合用於特徵萃取。橢圓形太扁,鄰近像素太少,無法充分萃取色彩變化。現在大家只找圓形blob。大家只做automatic scale selection,不做affine shape adaptation。

實作方式是利用Gaussian pyramid,採用特殊縮放倍率2¹⸍³ ≈ 1.26,仿照DoG filter的高速近似技巧,快速完成Harris corner detector與automatic scale selection。虛擬碼程式碼

Image Feature Description

特徵描述。圖片的鋩鋩角角,改寫成數值,以供辨識比對。

大家習慣一併考慮特徵萃取與特徵描述。

SIFT效果最好。ORB效果稍差但速度極快。

HOG   – histograms of oriented gradients
SIFT  – scale-invariant feature transform
GLOH  – gradient location and orientation histogram
SURF  – speeded up robust features
BRIEF – binary robust independent elementary features
ORB   – oriented FAST and rotated BRIEF
BRISK – binary robust invariant scalable keypoints
FREAK – fast retina keypoint
KAZE  – KAZE is a Japanese word that means wind
AKAZE – accelerated-KAZE

histograms of oriented gradients (HOG):一個特徵,調整成16×16正方形區域。一個區域256個像素實施投票,得到9個數字,RGB分開就是9×3個數字。一個區域的數字,依序串成一串,當作特徵描述。

一個像素進行投票,梯度長度是持票數,梯度角度(無視正反方向,僅0°到180°)決定投票箱編號。180°等分成9個箱子。一個像素投票給最鄰近的兩個箱子,與箱子中心的距離(角度差)作為投票比例。

scale-invariant feature transform (SIFT):一個特徵,調整成16×16正方形區域。一個區域切成4×4宮格。一個宮格所有像素實施投票,得到8個數字,所有宮格就是4×4×8=128個數字,RGB分開就是128×3個數字。一個區域的數字,依序串成一串,當作特徵描述。

一個像素進行投票,梯度長度是持票數,梯度角度決定投票箱編號。360°等分成8個箱子。一個像素投票給一個箱子。

我們希望圖片經過位移、縮放、旋轉,特徵描述仍然相同,達成相似不變。梯度具備位移不變。只看梯度相對長度,達成縮放不變:所有宮格所有投票箱得票數比例化(源自變種RootSIFT)。只看梯度相對角度,達成旋轉不變:一個宮格8個投票箱得票數循環移位,以得票數最多者為首。

gradient location and orientation histogram (GLOH):SIFT的改造版本。效果相仿但速度更慢。

speeded up robust features (SURF):SIFT的高速近似算法。效果稍差但速度較快。

binary robust independent elementary features (BRIEF):一個特徵,其鄰近像素形成41×41正方形區域。一個區域,任取兩個像素比較大小,得到0或1。重複128次,得到128個數字。一個區域的數字,依序串成一串,當作特徵描述。

oriented FAST and rotated BRIEF (ORB):追加了縮放不變與旋轉不變的FAST、追加了旋轉不變的BRIEF。區域中心、區域裡面每個像素位置的加權平均值(權重是像素值),兩者連線,旋轉對齊X軸,以便達成旋轉不變。

Image Feature Matching

特徵匹配。兩張圖片的鋩鋩角角,找到對應關係。

此處的圖片是臺北郵局臺北郵局,拍攝時間與取景地點略有差異。我為了讓匹配線條不交叉,特別縮小右圖尺寸。

最初大家只考慮特徵描述差異,後來大家也考慮特徵位置差距。

ANN最常使用。可能是因為參數數量最少,容易調參。

NN   – nearest neighbor
ANN  – approximate nearest neighbor
GMS  – grid-based motion statistics
CODE – coherence based decision boundaries

nearest neighbor (NN):每個圖一特徵,找到差異最小的圖二特徵。缺點是不保證一一對應。

approximate nearest neighbor (ANN):每個圖一特徵,找到差異接近最小的圖二特徵,只要誤差不超過ε倍即可。算得更快。知名函式庫FLANN

nearest neighbor distance ratio test:差異最小和差異次小的差距比值。比值越小越好。自訂臨界值。

cross check test:上述圖二特徵,找到描述差異接近最小的圖一特徵,判斷是否是相同圖一特徵。

grid-based motion statistics (GMS):一、事先篩選特徵,根據鄰近特徵數量多寡。二、每個特徵,以自己為中心,劃分九宮格,收錄鄰近特徵。以九宮格進行匹配。三、圖片縮放改成九宮格縮放,圖片旋轉改為九宮格輪轉,以達成縮放不變與旋轉不變。程式碼

coherence based decision boundaries (CODE):圖一到圖二的特徵移動方向,納入目標函數。其實就是homography alignment。程式碼

Image Feature Alignment

特徵對齊。兩張圖片的鋩鋩角角,找到變換函數。

最初大家只考慮特徵位置差距,後來大家也考慮特徵描述差異。

ICP最常使用。簡單粗暴效果好。

RANSAC – random sample consensus
PROSAC – progressive sample consensus
ICP    – iterative closest point
ICPIF  – ICP registration using invariant features

random sample consensus (RANSAC):隨機挑選多個特徵(數量固定,兩圖一樣多),隨機指定對應關係,找到最佳變換函數,求得位置差距總和。重複數個回合,取位置差距總和最小者。或者,重複無限回合,直到位置差距總和難以變小。

progressive sample consensus (PROSAC):袋子事先放入大量特徵,只從袋子隨機挑選多個特徵。袋子優先放入比較可靠的特徵。袋子每回合追加一個特徵。

iterative closest point (ICP):每個圖一特徵,找到位置差異接近最小的圖二特徵。找到最佳變換函數,套用最佳變換函數。反覆實施直到收斂。

ICP registration using invariant features (ICPIF):同時考慮特徵位置與特徵描述。

Image Alignment(Image Registration)

對齊對位。兩張圖片的所有像素,找到對應關係。

圖片實施image warping,盡量符合另一張圖片,找到變換函數。像素位置集體改變,像素數值盡量相符,求得改變方式。

大家習慣將所有像素改成少數特徵,節省時間,避開異物。問題變成image feature matching與image feature alignment,同時找到對應關係與變換函數。

Video Tracking🚧

簡介

影片的數學知識:Optical Flow

光流。影片的特定時刻(一張圖片),每個像素的運動速度。

此處的圖片源自這段影片。左圖是影片的特定時刻;想要求得光流,還需要下個時刻。中圖是光流,畫成箭矢圖,無法涵蓋每個像素。右圖也是光流,畫成濃度圖,只有速度大小,沒有速度方向。

一段影片,鏡頭和景物持續移動。每張圖片,像素位置持續改變。相鄰圖片,實施image alignment,得到像素運動方向暨距離。相鄰圖片時距已知,得到像素運動速度。相鄰圖片相差無幾,改為實施optical flow,降低計算時間。

下述兩式聯立,得到像素亮度梯度與像素速度的關係式。

一、假定像素亮度不變、位置會變。
二、假定像素直線運動(位移變換)。
⎧ (1) image alignment (x,y) -> (x+Δx,y+Δy)
⎪ I(x,y,t) = I(x+Δx,y+Δy,t+Δt)
⎨ 
⎪ (2) first-order Taylor approximation
⎪ I(x+Δx,y+Δy,t+Δt) ≈ I(x,y,t) + Iₓ(x,y,t)Δx
⎩                   + I₝(x,y,t)Δy + Iₜ(x,y,t)Δt

註:Δx是物理學displacement,不是數學finite difference!
Iₓ(x,y,t)Δx + I₝(x,y,t)Δy + Iₜ(x,y,t)Δt = 0       (2) - (1)

Iₓ(x,y,t)(Δx/Δt) + I₝(x,y,t)(Δy/Δt) + Iₜ(x,y,t) = 0   divide Δt

Iₓ(x,y,t)vₓ(t) + I₝(x,y,t)v₝(t) + Iₜ(x,y,t) = 0   velocity at t

Iₓ(x,y)vₓ + I₝(x,y)v₝ + Iₜ(x,y) = 0               skip t

[ Iₓ(x,y) I₝(x,y) ] ⎡ vₓ ⎤ = [ -Iₜ(x,y) ]         matrix form
                    ⎣ v₝ ⎦

利用此關係式,求得一個像素的運動速度。

一個像素形成一次方程式,鄰近像素們形成一次方程組。找到平方誤差最小的解,作為該像素的運動速度。

⎡ Iₓ(x₁,y₁) I₝(x₁,y₁) ⎤          ⎡ -Iₜ(x₁,y₁) ⎤
⎢ Iₓ(x₂,y₂) I₝(x₂,y₂) ⎥ ⎡ vₓ ⎤ = ⎢ -Iₜ(x₂,y₂) ⎥
⎢      :         :    ⎥ ⎣ v₝ ⎦   ⎢       :    ⎥
⎣ Iₓ(xₙ,yₙ) I₝(xₙ,yₙ) ⎦          ⎣ -Iₜ(xₙ,yₙ) ⎦
           A              x             b

⎡ ∑IₓIₓ ∑I₝Iₓ ⎤ ⎡ vₓ ⎤ = ⎡ -∑IₓIₜ ⎤
⎣ ∑IₓI₝ ∑I₝I₝ ⎦ ⎣ v₝ ⎦   ⎣ -∑I₝Iₜ ⎦
      AᵀA         x          Aᵀb

⎡ vₓ ⎤ = ⎡ ∑IₓIₓ ∑I₝Iₓ ⎤⁻¹ ⎡ -∑IₓIₜ ⎤
⎣ v₝ ⎦   ⎣ ∑IₓI₝ ∑I₝I₝ ⎦   ⎣ -∑I₝Iₜ ⎦
  x          (AᵀA)⁻¹           Aᵀb

x = (AᵀA)⁻¹Aᵀb   normal equation

一個像素:A獲得1個橫條,b獲得1個元素。
挑選至少2個像素,獲得至少2道式子,以便解出2個變數。

按理來說,鄰近像素越多,答案越精準。然而鄰近像素速度不見得完全一致。鄰近像素越多,答案越不精準。

解法是scaling method。縮放尺寸、疊加答案。Gaussian pyramid,從最小的圖片(最大的範圍)開始計算運動速度。先以上回合的運動速度來移動像素,再求得這回合的運動速度。運動速度逐漸細膩。稱作coarse-to-fine optical flow estimation。

Optical Flow Estimation

光流估計。影片特定時刻(一張圖片)的每個像素的運動速度。

Lucas–Kanade algorithmlinear least squares。「像素數值的差異」越小越好。即是上個小節的演算法。光流的發明人。

Horn–Schunck algorithmmultiobjective optimization。「像素數值的差異」與「速度向量的長度」越小越好。

EpicFlowmultiobjective optimization。「像素數值的差異」與「速度向量的長度」與「邊緣匹配的差異」越小越好。

Motion Estimation

運動估計。影片特定時刻(一張圖片)的物體運動速度。

物體運動跟像素運動不見得相同。例如球狀物體原地自轉。例如物體不動而光源移動。

影片物體運動估計,我們只能以像素運動當作物體運動。後面章節的相機物體運動估計,藉由多種視角,我們也許有更好的演算法,不過目前還沒有人研究這個主題。

Image Template Tracking

範本追蹤。給定圖片片段,找到移動軌跡。

Lucas–Kanade algorithm:圖片片段與第一張圖片做image alignment得到初始位置。相鄰兩兩圖片做optical flow estimation得到移動方向。銜接成移動軌跡。

Image Feature Tracking

特徵追蹤。根據圖片特徵,找到移動軌跡。

Kanade–Lucas–Tomasi algorithm:一、光流估計:全部像素改成部分特徵。位移變換改成仿射變換。公式解改成疊代法。詳細推導過程請見這篇講義。二、特徵偵測:Shi–Tomasi corner detector。二次動差矩陣總和的兩個特徵值的絕對值大小,原本用來判斷是否有角,現在還用來判斷運動速度是否有效。

Video Object Tracking

物體追蹤。指定圖片當中物體,找到移動軌跡。

進階應用有多物體追蹤群眾計數

tracking by alignment:template tracking (dense tracking)或feature tracking (sparse tracking),請見前面小節。

tracking by similarity:早期作法是online optimization,以intensity histogram或intensity-gradient historgram (HOG/SIFT)的相似程度作為函數,實施mean shift或particle filter。近期作法是deep learning,我沒有研究。

Camera Tracking (with object coordinates)🚧

簡介

Camera Imaging

成像。一臺相機,拍攝一個物體,得到一張相片。

pinhole imaging model:針孔成像。

物體發射光線,穿越針孔,抵達圖片。上下左右顛倒。

三個元件。

一、物體:幾何模型。請見本站文件Model。
二、相機:幾何變換。請見本站文件Transformation。
三、圖片:二維陣列。請見本站文件Image

原點是相機針孔,Z軸朝向圖片,右手座標系。

一條射線形成三個點座標。

object (x,y,z)
image  (u,v,w)
camera (0,0,0)

三個點座標的數學關係式。相似三角形對應邊長比例相等。

⎧ u = -(f/z) x
⎨ v = -(f/z) y
⎩ w = -(f/z) z = -f

數位相機組裝多種透鏡,不是針孔成像。然而無論組裝哪種透鏡,皆可視作針孔成像,唯有焦距f產生變化。針孔成像足夠堪用。

數位相機製造商自己知道詳細構造,不需要特地視作針孔成像。一般民眾無法得知詳細構造,不得已視作針孔成像。

【待補圖片】

perspective projection model:透視投影。

計算學家為了省略負號而搞出來的玩意兒。

圖片挪至針孔的另一側。不必上下左右顛倒。

⎧ u = (f/z) x
⎨ v = (f/z) y
⎩ w = (f/z) z = f

Tsai's camera model:區分相機內外。

三個元件、三種座標系。

三個步驟。

[相機外]一、物體座標,從物體座標系更換成相機座標系。
       座標軸的相似變換(旋轉、位移、縮放)

[相機內]二、物體座標變成圖片座標,根據相機成像機制。
       透視投影(非一次函數)
       扭曲(非一次函數)

     三、圖片座標,從相機座標系更換成圖片座標系。
       圖片的長寬(縮放)
       主點的索引值(位移)

三個步驟當中所有參數。分別稱作外在參數和內在參數。

[相機外]方向   orientation             R
     位置   position                tₓ t₝ t₞
     比例尺度 ratio scale             s

[相機內]焦距   focal length            f
     扭曲係數 distortion coefficient  k₁ k₂

     長寬比  aspect ratio            sᵤ sᵥ
     偏移量  offset                  oᵤ oᵥ

步驟一三可以寫成矩陣運算。步驟二無法寫成矩陣運算。

1. change coordinate system of object (object -> camera)

⎡ xꟲ ⎤     ⎡ R₁₁ R₁₂ R₁₃ tₓ ⎤ ⎡ xᴼ ⎤
⎢ yꟲ ⎥ = s ⎢ R₂₁ R₂₂ R₂₃ t₝ ⎥ ⎢ yᴼ ⎥     s can be eliminated
⎣ zꟲ ⎦     ⎣ R₃₁ R₃₂ R₃₃ t₞ ⎦ ⎢ zᴼ ⎥     due to perspective
                              ⎣ 1  ⎦     projection

2. transform object to image (in camera coordinate system)

⎰ uꟲ = (f/zꟲ) xꟲ                perspective projection
⎱ vꟲ = (f/zꟲ) yꟲ

⎰ úꟲ = (1 + k₁r² + k₂r⁴) uꟲ     distortion
⎱ v́ꟲ = (1 + k₁r² + k₂r⁴) vꟲ     where r² = (uꟲ)² + (vꟲ)²

3. change coordinate system of image (camera -> image)

⎡ uᴵ ⎤   ⎡ sᵤ  0 oᵤ ⎤ ⎡ úꟲ ⎤
⎣ vᴵ ⎦ = ⎣  0 sᵥ oᵥ ⎦ ⎢ v́ꟲ ⎥
                      ⎣ 1  ⎦

camera matrix model:省略扭曲,強行化作矩陣運算。

三式併成一式y = λ K [R|t] x,清爽暢快。

兩個矩陣分別稱作外在矩陣[R|t]和內在矩陣K。

⎡ uᴵ ⎤     ⎡ fᵤ  0 zᵤ ⎤ ⎡ R₁₁ R₁₂ R₁₃ tₓ ⎤ ⎡ xᴼ ⎤   where
⎢ vᴵ ⎥ = λ ⎢  0 fᵥ zᵥ ⎥ ⎢ R₂₁ R₂₂ R₂₃ t₝ ⎥ ⎢ yᴼ ⎥   fᵤ = f sᵤ
⎣ 1  ⎦     ⎣  0  0  1 ⎦ ⎣ R₃₁ R₃₂ R₃₃ t₞ ⎦ ⎢ zᴼ ⎥   fᵥ = f sᵥ
                                           ⎣ 1  ⎦   λ  = 1/zꟲ
image      intrinsic    extrinsic          object   zᵤ = zꟲ oᵤ
coord.     matrix       matrix             coord.   zᵥ = zꟲ oᵥ
y          K            [R|t]              x

矩陣運算的細節。

[外在矩陣]
一、旋轉R、位移t,讓物體自由調整擺放姿勢。不必煩惱原點水平線。
  縮放s,讓物體自由調整度量單位。不必煩惱公分英吋。
二、相似變換可以採用任意順序。原作者採用旋轉R、位移t、縮放s。
  如此一來,透視投影恰巧消掉縮放s。節省一個參數s。
三、位移不是線性函數,但是位移是一次函數。
  使用齊次座標(補1),得以寫成矩陣運算。
[內在矩陣]
一、內在矩陣K添加一個橫條,方便矩陣連乘。
  順帶一提,內在矩陣K恰是上三角矩陣。
二、透視投影、扭曲不是線性函數。
  無法寫成矩陣運算。
三、透視投影不是線性函數,但是透視投影是射影線性函數。
  使用射影座標(補1、等於=改成等價≡),得以佯裝成矩陣運算。
  本文不使用射影座標。(射影幾何不是這樣用的。)
  教科書使用射影座標。(齊次座標與射影座標外觀幾乎相同,經常導致混淆。)
[尺度參數λ]
一、透視投影硬是納入矩陣運算。
 甲、相機座標系的Z軸,垂直於圖片、穿過圖片中央(主點),其座標是焦距f。
   透視投影uꟲ = (f/zꟲ)xꟲ和vꟲ = (f/zꟲ)yꟲ。
 乙、乘以f是一次函數,挪到內在矩陣。
   併入sᵤ和sᵥ。
   新參數fᵤ = f sᵤ和fᵥ = f sᵥ。意義是pixel density。
 丙、除以zꟲ不是一次函數,挪到矩陣外面最左邊。
   新參數λ = 1/zꟲ。
   副作用是oᵤ和oᵥ必須隨之調整。
   新參數zᵤ = zꟲ oᵤ和zᵥ = zꟲ oᵥ。
二、λ挪到矩陣外面最左邊,得以佯裝成射影線性函數。
  然而λ = 1/zꟲ不是真正的尺度,畢竟矩陣內含z值。

Camera Lens Distortion

鏡頭扭曲。透鏡曲面不完美,透鏡組合不完美,投影位置偏移。尤其是魚眼鏡頭fisheye lens、全景相機omnidirectional camera。

Brown–Conrady model:多項式函數,變數是半徑平方。

radical distortion:
duᵣ = (k₁r² + k₂r⁴ + ...) u     where r = sqrt(u² + v²)
dvᵣ = (k₁r² + k₂r⁴ + ...) v     skip ꟲ for convenience

tangential distortion:
duₜ = p₁(r² + 2u²) + 2p₂uv
dvₜ = p₂(r² + 2v²) + 2p₁uv

distortion:
ú = u + duᵣ + duₜ
v́ = v + dvᵣ + dvₜ

Kannala–Brandt model:多項式函數,變數是入射角平方。

du = (1/r) (k₁θ + k₂θ³ + k₃θ⁵ + ...) u     where r = sqrt(u² + v²)
dv = (1/r) (k₁θ + k₂θ³ + k₃θ⁵ + ...) v           θ = tan⁻¹(r/f)

Camera Calibration / Camera Resectioning

校準、率定。一臺相機,求得相機成像參數。求得內在參數。

交會。一臺相機,找到相機相對於物體的方位。換句話說,找到物體相對於相機的位置。求得外在參數。

不知位置,無從校準。未經校準,不得位置。兩者相依。兩者必須一併求得。大家將兩者視作相同問題。

camera calibration
given ((xᴼ,yᴼ),(uᴵ,vᴵ)) find f k₁ k₂ sᵤ sᵥ oᵤ oᵥ

camera resectioning
given ((xᴼ,yᴼ),(uᴵ,vᴵ)) find R t

camera calibration = camera resectioning
given ((xᴼ,yᴼ),(uᴵ,vᴵ)) find R t f k₁ k₂ sᵤ sᵥ oᵤ oᵥ

成像視作函數求值,校準視作迴歸。一、測量物體。二、拍攝相片。三、找出物體與相片的對應點對(的座標)。擷取多組對應點對,提升迴歸準確度。四、求得相機成像機制所有參數。

數位相機製造商自己知道詳細構造,不需要特地求得成像參數。一般民眾無法得知詳細構造,不得已求得成像參數。

Tsai's algorithm:相機成像三個步驟分批處理。

已知第三步驟參數sᵤ sᵥ oᵤ oᵥ。先求第一步驟參數R tₓ t₝,再求第二步驟參數t₞ f k₁ k₂。

需要事先知道第三步驟參數sᵤ sᵥ oᵤ oᵥ。實務上是用圖片尺寸、感光元件尺寸求得sᵤ sᵥ oᵤ oᵥ。實務上需要通靈問出尺寸

sᵤ = image_width / sensor_width
sᵥ = image_height / sensor_height
oᵤ = image_width / 2
oᵥ = image_height / 2

未命名演算法:傳統工藝菁華集成,內容源自多篇論文。

一、物體與圖片的對應點對:解法有自訂標記、西洋棋盤。

一甲、image feature detection。物體事先塗上數個特殊點,方便找到對應點對。特殊點事先測量座標,方便計算成像參數。

一乙、chessboard detection。三個西洋棋盤拼成凹陷立方體,方格四角當作特殊點。不必事先塗上特殊點,不必事先測量座標,只需實施直線偵測,只需實施線段相交。

二、內在外在參數:nonlinear least squares。疊代法有block coordinate descent、BFGS Algorithm、Levenberg–Marquardt algorithm。大家採用後者。

   min   sum ‖yᵢ - F(xᵢ; R,t,λ,f,⋯)‖²
R,t,λ,f,⋯ i

三、疊代法初始值:差不多就行了。因此採用內在外在矩陣。

三甲、內在外在矩陣乘積:homogeneous linear least squares

λ挪到矩陣外面最左邊,硬是佯裝成homography alignment。消掉λ之手法即是direct linear transformation。消掉λ所形成的兩道一次方程式,作者稱作collinearity equation。答案是AᵀA的最小的特徵值的特徵向量。

因為省略扭曲,所以只能用來估計初始值。

此步驟源自Abdel-Aziz的這篇論文Sutherland的這篇論文

⎡ uᴵ ⎤     ⎡ H₁₁ H₁₂ H₁₃ H₁₄ ⎤ ⎡ xᴼ ⎤
⎢ vᴵ ⎥ = λ ⎢ H₂₁ H₂₂ H₂₃ H₂₄ ⎥ ⎢ yᴼ ⎥
⎣ 1  ⎦     ⎣ H₃₁ H₃₂ H₃₃ H₃₄ ⎦ ⎢ zᴼ ⎥
                               ⎣ 1  ⎦
image      product             object
coord.     of two matrices     coord.
y          H = K [R|t]         x

⎧ λ (H₁₁x + H₁₂y + H₁₃z + H₁₄) = u     skip ᴵ ᴼ for convenience
⎨ λ (H₂₁x + H₂₂y + H₂₃z + H₂₄) = v
⎩ λ (H₃₁x + H₃₂y + H₃₃z + H₃₄) = 1

⎰ (H₁₁x + H₁₂y + H₁₃z + H₁₄) - u (H₃₁x + H₃₂y + H₃₃z + H₃₄) = 0
⎱ (H₂₁x + H₂₂y + H₂₃z + H₂₄) - v (H₃₁x + H₃₂y + H₃₃z + H₃₄) = 0

                                   ⎡ H₁₁ ⎤
⎡ x y z 1 0 0 0 0 -ux -uy -uz -u ⎤ ⎢ H₁₂ ⎥ = ⎡ 0 ⎤
⎣ 0 0 0 0 x y z 1 -vx -vy -vz -v ⎦ ⎢  :  ⎥   ⎣ 0 ⎦
                                   ⎣ H₃₄ ⎦
                A                     x        0

solve Ax = 0   --->   min ‖Ax‖² = 0 subject to ‖x‖² = 1

x = smallest eigenvector of AᵀA
  = smallest singular vector of AᵀA
    (since AᵀA is positive semidefinite)

一組對應點對:A獲得2個橫條,零向量獲得2個元素。
挑選至少6組對應點對,獲得至少12道式子,以便解出12個變數。

三乙、內在外在矩陣:RQ decomposition。

因為K₁₂不見得是零(歪斜),所以只能用來估計初始值。

此步驟源自Hartley的這篇論文。多重眼球書的作者。

A = 𝑅𝑄             RQ decomposition
A = 𝑄′𝑅′           QR decomposition

𝑅 and 𝑅′ are upper triangular matrices
𝑄 and 𝑄′ are orthonormal matrices
get “RQ decomposition of A” by “QR decomposition of A⁻¹”

A = 𝑅𝑄             RQ decomposition of A   (unknown)
A⁻¹ = 𝑄″𝑅″         QR decomposition of A⁻¹ (known)
A = 𝑅″⁻¹𝑄″⁻¹       inverse formula

⎧ 𝑅 = 𝑅″⁻¹         inverse of upper triangular matrix
⎨                  is also an upper triangular matrix
⎪ 𝑄 = 𝑄″⁻¹ = 𝑄″ᵀ   inverse of orthonormal matrix
⎩                  is also an orthonormal matrix
H = K [R|t]

⎡ |  |  |  |  ⎤     ⎡ |  |  |  | ⎤
⎢ h₁ h₂ h₃ h₄ ⎥ = K ⎢ r₁ r₂ r₃ t ⎥
⎣ |  |  |  |  ⎦     ⎣ |  |  |  | ⎦

⎡ |  |  |  ⎤
⎢ h₁ h₂ h₃ ⎥ = K R = 𝑅 𝑄
⎣ |  |  |  ⎦

⎧ K = 𝑅
⎨ R = 𝑄
⎩ t = K⁻¹h₄

resolve name collision:
R is rotation matrix (orthonormal matrix)
𝑅 is from RQ decomposition (upper triangular matrix)

三丙、尺度參數:maximum likelihood estimation。考慮每組對應點對之每個維度的尺度,求期望值(平均數)。

通常只挑選一組對應點對,節省時間。

y = λ K [R|t] x               camera imaging

y = λ H x                     product of two matrices

⎡ ŭ ⎤     ⎡ u ⎤
⎢ v̆ ⎥ = λ ⎢ v ⎥               vector form
⎣ 1 ⎦     ⎣ w ⎦

λ = (ŭ/u + v̆/v + 1/w) / 3     expected value
⎧ y₁ = λ₁ H x₁                camera imaging of N points
⎨  :      :
⎩ yɴ = λɴ H xɴ

λ = (λ₁ + ... + λɴ) / N       expected value (slow)

(sum yᵢ) = λ H (sum xᵢ)       additivity

⎡ ŭ ⎤     ⎡ u ⎤
⎢ v̆ ⎥ = λ ⎢ v ⎥               vector form
⎣ 1 ⎦     ⎣ w ⎦

λ = (ŭ/u + v̆/v + 1/w) / 3     expected value (fast)

三丁、扭曲係數:linear least squares。公式解有normal equation、QR decomposition、SVD。作者採用前者。

此步驟源自Zhang的這篇論文。下一個介紹的演算法。

Brown–Conrady model
⎰ úꟲ = (1 + k₁r² + k₂r⁴) uꟲ     where r² = (uꟲ)² + (vꟲ)²
⎱ v́ꟲ = (1 + k₁r² + k₂r⁴) vꟲ

difference
⎰ (úꟲ - uꟲ) = (k₁r² + k₂r⁴) uꟲ
⎱ (v́ꟲ - vꟲ) = (k₁r² + k₂r⁴) vꟲ

let uᴵ vᴵ be the estimation from camera imaging y = λ K [R|t] x
let ŭᴵ v̆ᴵ be the observation from image feature detection
⎰ (1/sᵤ) (ŭᴵ - uᴵ) = (k₁r² + k₂r⁴) (1/sᵤ) (uᴵ - oᵤ)
⎱ (1/sᵥ) (v̆ᴵ - vᴵ) = (k₁r² + k₂r⁴) (1/sᵥ) (vᴵ - oᵥ)
where r² = ((1/sᵤ) (uᴵ - oᵤ))² + ((1/sᵥ) (vᴵ - oᵥ))²

linear least squares
⎡ (uᴵ-oᵤ)r² (uᴵ-oᵤ)r⁴ ⎤ ⎡ k₁ ⎤ = ⎡ ŭᴵ-uᴵ ⎤
⎣ (vᴵ-oᵥ)r² (vᴵ-oᵥ)r⁴ ⎦ ⎣ k₂ ⎦   ⎣ v̆ᴵ-vᴵ ⎦
           A              x          b

一組對應點對:A獲得2個橫條,b獲得2個元素。
挑選至少1組對應點對,獲得至少2道式子,以便解出2個變數。

四、疊代法過程:不採用內在外在矩陣。

四甲、旋轉:可以表示成旋轉矩陣(9個變數)、自旋四元數(4個變數)、自旋四元數對數(3個變數)。大家採用後者。

精簡變數以降低疊代法計算時間。精簡變數以避免約束條件。3個變數恰是旋轉自由度,不能再少了。

some formulas for quaternion
q = a + bi + cj + dk
exp(q) = exp(a) [ cos(‖v‖) + (v/‖v‖) sin(‖v‖) ]
ln(q) = ln(‖q‖) + (v/‖v‖) cos⁻¹(a/‖q‖)

let real part be a
let imagine part be v = bi + cj + dk
quaternion of spin
q = exp(u ½θ)
  = cos(½θ) + u sin(½θ)
  = cos(½θ) + uₓ sin(½θ) i + u₝ sin(½θ) j + u₞ sin(½θ) k

logarithm of quaternion of spin (real part is zero!)
p = ln(q) = u ½θ

let spin axis be unit vector (uₓ,u₝,u₞)   (uₓ²+u₝²+u₞² = 1)
let spin angle be θ
let quaternion of spin axis be u = uₓi + u₝j + u₞k
logarithm of quaternion of spin
p = wₓi + w₝j + w₞k

quaternion of spin
exp(p) = exp(0) [ cos(‖w‖) + (w/‖w‖) sin(‖w‖) ]
       = cos(½θ) + u sin(½θ)

let real part be 0
let imagine part be w = wₓi + w₝j + w₞k
thus w/‖w‖ is spin axis
thus 2‖w‖ is spin angle
cross product matrix
       ⎡   0 -w₞  w₝ ⎤
[w]ₓ = ⎢  w₞   0 -wₓ ⎥
       ⎣ -w₝  wₓ   0 ⎦

rotation matrix (Rodrigues' formula)
        sin(2‖w‖)        1 - cos(2‖w‖)
R = I + ————————— [w]ₓ + ————————————— [w]ₓ²
           ‖w‖                ‖w‖²

四乙、焦距:初始值算的是fᵤ fᵥ。我們想算的是f sᵤ sᵥ。

需要事先知道sᵤ sᵥ,才能將fᵤ fᵥ分解為f sᵤ sᵥ。實務上是用圖片尺寸、感光元件尺寸求得sᵤ sᵥ。實務上需要通靈問出尺寸

sᵤ = image_width / sensor_width
sᵥ = image_height / sensor_height
f = fᵤ / sᵤ = fᵥ / sᵥ

然而大家沒有分解。

疊代法過程,繼續使用fᵤ fᵥ。疊代法結束,才分解成f sᵤ sᵥ。

由於省略了扭曲,這是錯誤的作法!但是大家都用得很開心,一片和諧。OpenCVMATLAB皆是如此。

四丙、目標函數的梯度(Jacobian matrix):硬著頭皮手工推導吧。一片和諧。

Zhang's algorithm:改用二維物體。主流演算法。

一、物體:置於多種方位,拍攝多張相片,避免共面導致多解。

「一個三維物體」改成「多個二維物體」。
「多個二維物體」改成「同一個二維物體置於多種方位」。
「聯合拍攝相片」改成「獨自拍攝相片」。

二、物體與圖片的對應點對:chessboard detection

三、疊代法:同前。針對多張圖片,修改計算流程。

N張圖片,每張圖片抓M組對應點對,加總所有誤差。

           N   M
   min    sum sum ‖yₙᵢ - Fₙ(xₙᵢ; R,t,λ,f,⋯)‖²
R,t,λ,f,⋯ n=1 i=1

四、疊代法初始值:針對二維物體,稍微修改計算流程。

四甲、內在外在矩陣乘積:homogeneous linear least squares。同前。每張圖片的外在矩陣都不同。

四乙、內在矩陣倒數內積:homogeneous linear least squares。旋轉矩陣恰是正規正交矩陣,以此製造等式。

作者並未活用「向量長度等於一」的性質。作者令B = λK⁻ᵀK⁻¹,憑空多出λ,要嘛也是λ²。我認為推導過程出現紕繆!但是大家都用得很開心,一片和諧。

⎡ uᴵ ⎤     ⎡ fᵤ  0 zᵤ ⎤ ⎡ |  |  |  | ⎤ ⎡ xᴼ ⎤
⎢ vᴵ ⎥ = λ ⎢  0 fᵥ zᵥ ⎥ ⎢ r₁ r₂ r₃ t ⎥ ⎢ yᴼ ⎥   2D object
⎣ 1  ⎦     ⎣  0  0  1 ⎦ ⎣ |  |  |  | ⎦ ⎢ 0ᴼ ⎥   z = 0
                                       ⎣ 1  ⎦
image      intrinsic    extrinsic      object
coord.     matrix       matrix         coord.
y          K            [R|t]          x

⎡ uᴵ ⎤     ⎡ fᵤ  0 zᵤ ⎤ ⎡ |  |  | ⎤ ⎡ xᴼ ⎤
⎢ vᴵ ⎥ = λ ⎢  0 fᵥ zᵥ ⎥ ⎢ r₁ r₂ t ⎥ ⎢ yᴼ ⎥   2D object
⎣ 1  ⎦     ⎣  0  0  1 ⎦ ⎣ |  |  | ⎦ ⎣ 1  ⎦   z = 0

        ⎡ |  |  |  ⎤     ⎡ |  |  | ⎤
let H = ⎢ h₁ h₂ h₃ ⎥ = K ⎢ r₁ r₂ t ⎥
        ⎣ |  |  |  ⎦     ⎣ |  |  | ⎦
⎧ r₁ = K⁻¹h₁
⎨ r₂ = K⁻¹h₂
⎪ r₁∙r₂ = 0             orthonormal matrix: perpendicularity
⎩ r₁∙r₁ = r₂∙r₂ = 1     orthonormal matrix: unit length

let r₁∙r₁ = r₁ᵀr₁
⎰ h₁ᵀK⁻ᵀK⁻¹h₂ = 0
⎱ h₁ᵀK⁻ᵀK⁻¹h₁ - h₂ᵀK⁻ᵀK⁻¹h₂ = 0

let B = K⁻ᵀK⁻¹          symmetric matrix
⎰ h₁ᵀBh₂ = 0
⎱ h₁ᵀBh₁ - h₂ᵀBh₂ = 0

let v₁ = [ H₁₁H₁₂ H₂₁H₂₂ H₃₁H₃₂ H₁₁H₂₂ H₁₁H₃₂ H₂₁H₃₂ ]
       + [      0      0      0 H₂₁H₁₂ H₃₁H₁₂ H₃₁H₂₂ ]

let v₂ = [ H₁₁² H₂₁² H₃₁² 2H₁₁H₂₁ 2H₁₁H₃₁ 2H₂₁H₃₁ ]
       - [ H₁₂² H₂₂² H₃₂² 2H₁₂H₂₂ 2H₁₂H₃₂ 2H₂₂H₃₂ ]

             ⎡ B₁₁ ⎤
             ⎢ B₂₂ ⎥
⎡ —— v₁ —— ⎤ ⎢ B₃₃ ⎥ = ⎡ 0 ⎤
⎣ —— v₂ —— ⎦ ⎢ B₁₂ ⎥   ⎣ 0 ⎦
             ⎢ B₁₃ ⎥
             ⎣ B₂₃ ⎦
     A          x        0

solve Ax = 0   --->   min ‖Ax‖² = 0 subject to ‖x‖² = 1

x = smallest eigenvector of AᵀA
  = smallest singular vector of AᵀA
    (since AᵀA is positive semidefinite)

一張圖片(一種外在矩陣):A獲得2個橫條,零向量獲得2個元素。
挑選至少3張圖片(3種外在矩陣),獲得至少6道式子,以便解出6個變數。

四丙、內在矩陣:conjugate decomposition。解法有特徵分解、Cholesky分解。然後再求反矩陣。作者手工推導公式解。

因為K₁₂不見得是零(歪斜),所以只能用來估計初始值。

⎧ zᵥ = (B₁₂ B₁₃ - B₁₁ B₂₃) / (B₁₁B₂₂ - B₁₂²)
⎪ λ  = B₃₃ - [B₁₃² + zᵥ (B₁₂ B₁₃ - B₁₁ B₂₃)] / B₁₁
⎨ fᵤ = sqrt(λ / B₁₁)
⎪ fᵥ = sqrt(λ B₁₁ / (B₁₁ B₂₂ - B₁₂²))
⎪ zᵤ = - B₁₃ fᵤ² / λ
⎩ K₁₂ = - B₁₂ fᵤ² fᵥ / λ

四丁、外在矩陣:linear equation。H = K [R|t]移項即得。

⎧ λ  = 1/‖K⁻¹h₁‖ = 1/‖K⁻¹h₂‖
⎪ r₁ = λ K⁻¹ h₁
⎨ r₂ = λ K⁻¹ h₂
⎪ r₃ = r₁ × r₂
⎩ t  = λ K⁻¹ h₃

四戊、扭曲係數:linear least squares。同前。N張圖片,每張圖片抓M組對應點對,形成2NM個橫條。

五、疊代法過程:作者沒有提及,一片和諧。

Camera Pose Estimation

姿態估計。已知物體座標,找到相機相對於物體的方位。

換句話說,找到物體相對於相機的位置。求得外在參數。

camera calibration
given ((xᴼ,yᴼ),(uᴵ,vᴵ)) find R t f k₁ k₂ sᵤ sᵥ oᵤ oᵥ

camera pose estimation
given ((xᴼ,yᴼ),(uᴵ,vᴵ)) find R t

camera pose estimation (with calibrated camera) ✔
given ((xᴼ,yᴼ),(uᴵ,vᴵ)) f k₁ k₂ sᵤ sᵥ oᵤ oᵥ find R t

如果事先知道內在參數,那麼問題可以簡化、演算法可以簡化:

一、已知圖片座標、物體座標、對應點對。圖片座標更換成相機座標系,進行扭曲反運算。問題化作:物體座標如何更換成相機座標系(相似變換),使得物體座標經過透視投影符合圖片座標。

二、透視投影恰巧消掉縮放s。相似變換R t s簡化成剛體變換R t。問題化作:一群點,經過剛體變換、透視投影,對齊另外一群點。已知兩群點對應方式、透視投影方式。找到剛體變換。

三、已知內在參數,求得外在參數。未知變數降為6個(旋轉3個、位移3個)。挑選3組對應點對,獲得6道式子,即得唯一解。挑選多組對應點對,可得最小平方解。

計算學家將倆問題稱作perspective-3-point problem (P3P)和perspective-n-point problem (PnP)。相似問題如下。

orthogonal procrustes problem:找到剛體變換/相似變換。
演算法是Kabsch–Umeyama algorithm。

homography alignment:找到射影變換。
演算法是direct linear transformation。

perspective-n-point problem:找到剛體變換,然後透視投影。
演算法是iterative closest point。

Perspective-3-point Problem (P3P)

三個點,經過剛體變換、透視投影,對齊另外三個點。已知兩群點對應方式、透視投影方式。找到剛體變換。

P3P已經發展許多演算法。發展歷程大致是:矩陣運算、多項式運算、手工推導公式解。計算時間越來越短。

direct linear transformation:矩陣運算。沿用先前手法。

Grunert's algorithm:多項式運算。考慮物體座標系。兩個物體座標(已知)、焦點座標(未知),形成三角形。焦點座標到兩個物體座標的射線距離(未知)與射線夾角(已知)(從相機座標系的圖片座標求得射線夾角),兩個物體座標的距離(已知),總共三邊一角,形成餘弦定理。3個物體座標,3個三角形,3個餘弦定理,以此製造等式。重新整理為4次多項式方程式,變數是射線距離比值。解出射線距離,求焦點座標,求射線向量,求剛體變換。4次多項式方程式有4組解,其中一個是符合剛體變換的正解。

lambda twist:手工推導公式解。我沒有研究。

Perspective-n-point Problem (PnP)

一群點,經過剛體變換、透視投影,對齊另外一群點。已知兩群點對應方式、透視投影方式。請找到剛體變換。

PnP目前沒人找到公式解。大家採用最佳化。

Lu–Hager–Mjolsness algorithmiterative closest point。找到最佳剛體變換。為了應付透視投影,令誤差是點到射線的最短距離。

EPnPhomogeneous linear least squares。採用重心座標,捨棄連線距離。重心座標經過透視投影,座標仍舊不變,以此製造等式。追加約束條件「控制點的相對距離仍舊不變」,以得到唯一解。答案退化成0/1/2/3維,分開討論4種情況。詳情我沒有研究。

SQPnPsequential quadratic programming。我沒有研究。

Camera Motion Estimation

運動估計。相機運動,形成軌跡。拍攝多張相片,源自各種方位。求出每張相片的相機方位。

每張相片各自進行camera pose estimation。目前大家並未考慮相片之間的關係。既然有個座標明確的物體,與其參考鄰近相片,不如參考物體。

Camera Tracking (without object coordinates)🚧

簡介

Two-view Imaging

雙視角成像。兩臺相機,拍攝一個物體,得到兩張相片。

triangulation:三角測量。

一個量測目標,兩個測量地點,形成三角形。

一、已知測量地點距離、量測夾角,可得目標距離。

國中數學,ASA全等,兩角夾一邊,決定三角形造型。

高中數學,正弦定理,已知兩角夾一邊,可得三角三邊。

二、又知測量地點座標、量測仰角,可得目標座標。

高中數學,球座標系,已知半徑、經度、緯度,可得直角座標。

高中數學,射線前進,已知起點、方向、距離,可得終點。

two-view imaging model:相機成像暨三角測量。

一臺相機無法決定物體座標。不知深淺。

兩臺相機得以決定物體座標。三角測量。

一個物體,兩臺相機,形成三角形。

有三種描述方式,可以得到物體座標。

一、已知相機壹座標c₁、三個向量(相機壹到相機貳c₁c₂—→、相機壹到圖片壹c₁y₁—→、相機貳到圖片貳c₂y₂—→),可得物體座標。

二、已知相機壹射線(c₁, c₁y₁—→)、相機貳射線(c₂, c₂y₂—→),可得物體座標。

三、已知相機壹座標系之圖片壹座標y₁ꟲ¹=(u₁ꟲ¹,v₁ꟲ¹,f₁)、相機貳座標系之圖片貳座標y₂ꟲ²=(u₂ꟲ²,v₂ꟲ²,f₂)、相機貳座標系到相機壹座標系的轉換公式s[R|t],可得相機壹座標系之物體座標。因為只需射線方向,所以可以消掉縮放s。

epipolar constraint:省略扭曲,強行化作矩陣運算。省略三角形造型,強行製造等式。

一、製造等式。(y₁ꟲ¹)ᵀ E (y₂ꟲ²) = 0

使用兩角一邊製造等式,等式包含三角函數,耗費計算時間。

改用三邊共面與相對方位製造等式,等式不含三角函數,節省計算時間。然而等式不夠精確嚴謹,無法決定三角形造型。

c₁y₁———→ ∙ (c₁c₂———→ × c₂y₂———→) = 0      coplanarity
y₁ꟲ¹ ∙ (t × Ry₂ꟲ²) = 0        relative pose
(y₁ꟲ¹)ᵀ [t]ₓ R (y₂ꟲ²) = 0     cross product matrix
(y₁ꟲ¹)ᵀ E (y₂ꟲ²) = 0          essential matrix
y₁ᵀ E y₂ = 0                  skip ꟲ¹ ꟲ²

三個向量共面,以此製造等式。兩個向量的叉積,得到三角形法向量。剩下那個向量與三角形法向量的點積,得到零。

採用相機壹座標系。c₁和y₁仍是c₁和y₁。c₁是原點(0,0,0),c₁y₁—→是y₁。c₂和y₂從相機貳座標系更換成相機壹座標系是t和Ry₂ꟲ²+t。c₂是t,c₂y₂—→是兩者相減Ry₂ꟲ²。

二、essential matrix。E = [t]ₓ R。

等式中間兩個矩陣的相乘結果,稱作基本矩陣E。這個名稱不知如何吐槽。總之,這個矩陣用來縮短數學式子,清爽暢快。

三、等式更換座標系。(y₁ᴵ¹)ᵀ F (y₂ᴵ²) = 0。

圖片座標從相機座標系更換成圖片座標系。省略扭曲。

⎧ y₁ᴵ¹ = K₁ d₁(y₁ꟲ¹)               1st camera imaging
⎨ y₂ᴵ² = K₂ d₂(y₂ꟲ²)               2nd camera imaging
⎩ (y₁ꟲ¹)ᵀ E (y₂ꟲ²) = 0             essential matrix

⎧ y₁ᴵ¹ = K₁ y₁ꟲ¹                   no distortion d₁()
⎨ y₂ᴵ² = K₂ y₂ꟲ²                   no distortion d₂()
⎩ (y₁ꟲ¹)ᵀ E (y₂ꟲ²) = 0             essential matrix

(K₁⁻¹ y₁ᴵ¹)ᵀ E (K₂⁻¹ y₂ᴵ²) = 0
(y₁ᴵ¹)ᵀ (K₁⁻ᵀ E K₂⁻¹) (y₂ᴵ²) = 0
(y₁ᴵ¹)ᵀ F (y₂ᴵ²) = 0               fundamental matrix
y₁ᵀ F y₂ = 0                       skip ᴵ¹ ᴵ²

四、fundamental matrix。F = K₁⁻ᵀ E K₂⁻¹。

等式中間一串矩陣的相乘結果,稱作基本矩陣F。這個名稱不知如何吐槽。總之,這個矩陣用來縮短數學式子,清爽暢快。

五、已知R t,求得E。E = [t]ₓ R。

按照定義。

六、已知E,求得R t。[R|t] = [UW±1Vᵀ|±u₃]。

極分解,然後湊出叉積矩陣。四組候選解,兩臺相機分別面對物體或背對物體。其中一組是正解,兩臺相機皆面對物體。

E = UΣVᵀ = (UΣWᵀUᵀ)(UWVᵀ) = [t]ₓR

    ⎡ |  |  |  ⎤ ⎡ 1 0 0 ⎤ ⎡ |  |  |  ⎤ᵀ
E = ⎢ u₁ u₂ u₃ ⎥ ⎢ 0 1 0 ⎥ ⎢ v₁ v₂ v₃ ⎥      full SVD
    ⎣ |  |  |  ⎦ ⎣ 0 0 0 ⎦ ⎣ |  |  |  ⎦

    ⎡ 0 -1 0 ⎤
W = ⎢ 1  0 0 ⎥     for making cross product matrix
    ⎣ 0  0 1 ⎦     (orthonormal matrix W⁻¹ = Wᵀ)

⎰ R = UWVᵀ or UWᵀVᵀ
⎱ t = u₃ or -u₃

Camera Auto-calibration(Camera Self-calibration)

自校準。一臺相機,求得相機成像參數。不知物體座標。

epipolar constraint,改寫成fundamental matrix,不省略扭曲。兩張圖片擷取多組對應點對,實施nonlinear least squares。

就這樣。順手補充兩個古代過時作法。可以用來估計初始值。

Hartley's algorithm:三張圖片。奇異值分解,得到相機成像參數。推導過程隱含了Kruppa's equation,只不過是將上個章節的推導過程換句話說而已。

Self-Calibration of Stationary Cameras
原始論文。

Kruppa's Equations Derived from the Fundamental Matrix
解釋Kruppa's Equation。

Tomasi–Kanade factorization method:單張圖片。奇異值分解,得到相機成像參數,甚至得到物體座標。然而單張圖片不可能正確重建物體,例如視覺錯覺

Shape and Motion from Image Streams under Orthography: A Factorization Method
相機成像模型,透視投影硬是改成仿射變換,形成一次函數,以便使用奇異值分解。

Factorization Methods for Projective Structure and Motion
相機成像模型,仿射變換硬是調整成透視投影。

PowerFactorization: 3D Reconstruction with Missing or Uncertain Data
走火入魔。

Camera Relative Pose Estimation

相對姿態估計。找到相機貳相對於相機壹的方位。

上策求得方位cy c。中策求得方位變換R t、再求cy c。下策求得基本矩陣E、再求R t、再求cy c。

目前大家一律採用下策。原因也許是(y₁ꟲ¹)ᵀ E (y₂ꟲ²) = 0清爽暢快。然而,epipolar constraint無法決定三角形造型,essential matrix estimation只能勉強得到一種造型。

essential matrix estimation是邪教。難道沒有其他方法了嗎?

Essential Matrix Estimation

已知圖片座標(相機座標系),求得基本矩陣。

事前準備:
兩臺相機,拍攝一個物體,得到兩張相片。
image feature description得到特徵描述。
image feature matching得到對應點對。
camera calibration或camera auto-calibration得到內在參數。
圖片座標從圖片座標系更換成相機座標系。
進行扭曲反運算。

6-point algorithm:homogeneous linear least squares。normalized direct linear transformation。epipolar constraint。

                                      ⎡ E₁₁ ⎤
                                      ⎢ E₁₁ ⎥
[ u₁u₂ u₁v₂ u₁ v₁u₂ v₁v₂ v₁ u₂ v₂ 1 ] ⎢ E₁₃ ⎥ = 0
                                      ⎢  :  ⎥
                                      ⎣ E₃₃ ⎦

Nistér's 5-point algorithm:polynomial equation。5組對應點對,製造10次多項式方程式。高速算法

Fundamental Matrix Estimation

已知圖片座標(圖片座標系),求得基本矩陣。

事前準備:
兩臺相機,拍攝一個物體,得到兩張相片。
image feature description得到特徵描述。
image feature matching得到對應點對。
優點:可以應付各種來源的相片。
缺點:因為省略扭曲,所以只能用來估計初始值。

Hartley's 8-point algorithm:homogeneous linear least squares。normalized direct linear transformation。SVD with rank(F) = 2。穩健算法

7-point algorithm:polynomial equation。7組對應點對,製造14次多項式方程式。

Essential Matrix Estimation備註

epipolar constraint無法決定三角形造型。

實務上,大家多取幾組對應點對,盡量改善答案。宛如通靈。

理論上,大家取最少組對應點對,剛好形成答案。彷彿邪教。

方才介紹的演算法都是邪教。

Fundamental Matrix Estimation備註

data normalization可以減少浮點數誤差。

至於此處data normalization發揮神效,原因是扭曲。

大家省略扭曲,求得錯誤答案。data normalization重新定位扭曲中心。遠離扭曲中心的點,變得靠近扭曲中心。這些點集合更加擬合透視投影,錯誤答案更加準確。

要是不省略扭曲,求得正確答案。我相信,無論是否做了data normalization,正確答案幾乎不受影響。

Camera Motion Estimation

運動估計。相機運動,形成軌跡。拍攝多張相片,源自各種方位。求出每張相片的相機方位,將相機方位整合到世界座標系。大家習慣將第一個相機座標系指定為世界座標系。

一、本回合兩張圖片,
  找到對應點對,求得基本矩陣,求得R⁽ⁿ⁾ t⁽ⁿ⁾。(第n個相機座標系)
二、本回合兩張圖片、上回合兩張圖片,
  計算對應點對的距離比例,以此微調t的長度。
  ‖t⁽ⁿ⁾‖ / ‖t⁽ⁿ⁻¹⁾‖ = ‖y₁⁽ⁿ⁾ - y₂⁽ⁿ⁾‖ / ‖y₁⁽ⁿ⁻¹⁾ - y₂⁽ⁿ⁻¹⁾‖
三、累計R t,求得cy c。(化作第一個相機座標系,即是世界座標系。)

epipolar constraint無法決定三角形造型。essential matrix estimation只能勉強得到一種造型。必須重新調整t,才能得到合適的相機方位。

對應點對的距離比例,毫無根據。

essential matrix estimation是邪教。難道沒有其他方法了嗎?

Scene Tracking🚧

簡介

Object Pose Estimation

姿態估計。找到物體中心位置(座標)與正面方向(向量)。

有人在標題開頭添加6DoF或者6D,讓人一聽就知道是位置(x,y,z)與方向(i,j,k)六個變數。

正常來說,先有形狀估計,才有姿態估計。不過最近流行deep learning,大筆一揮直接得到結果,先後順序無關緊要了。

https://github.com/ZhongqunZHANG/awesome-6d-object
https://github.com/YoungXIAO13/ObjectPoseEstimationSummary

Object Shape Estimation

形狀估計。找到物體表面位置(點雲)與表面方向(法向量)。

進行形狀估計,細分為許多種工具。

目前最精準的工具是光達。常見的測量儀器是雷射測距儀實際範例

本文介紹的工具是相機,可以用在測距儀失效的場合,例如雨天行駛的自駕車。實際範例

instrument | principle
----------------------------
lidar      | time of flight
camera     | light intensity

使用相機進行形狀估計,又細分為許多種作法。

目前最精準的作法是shape from distortion。常見的測量儀器是結構光掃描儀實際範例

本文介紹的作法是shape from correspondence,可以用在掃描儀失效的場合,例如大型建築。實際範例

approach                  | popular technique     | key
------------------------------------------------------------
shape from distortion     | structured light      | 表面彎曲
shape from silhouette     | visual hull           | 投影輪廓
shape from shading        | photometric stereo    | 顏色漸層
shape from polarization   |                       | 光線路徑
shape from focus          |                       | 對焦距離
shape from correspondence | structure from motion | 對應點對

structured light:物體照射特殊造型的光源,例如直線條碼。根據直線條碼的彎曲程度,推導物體表面走向。

visual hull:一張相片的物體範圍,形成角錐。每張相片的角椎的交集,形成物體的大致輪廓:保留凸面、保留鞍面、填平凹面。注意到,此輪廓不是凸多面體、不是凸包。

photometric stereo:相機固定方位,光源改變方位,拍攝多張相片。甚至光源固定方位,只拍攝一張相片。根據光線反射機制、物體色彩變化,推導物體表面法向量。教學影片教學文章

shape from polarization:承上,使用濾鏡篩選光線方向,讓答案更加精準。

shape from focus:窮舉焦距,拍攝多張相片。找出對焦準確、畫面清晰之處,根據焦距大小、成像機制,推導該處座標。

structure from motion:拍攝多張相片,找到對應點對。根據三角測量原理、相機成像機制,推導對應點對該處座標。

Scene Reconstruction(Structure From Motion)

重建。拍攝多張相片/一段影片,建立物體模型。

流程大致如下:

structure from motion
1. image feature extraction         圖片特徵萃取
2. image feature description        圖片特徵描述
3. image feature matching           圖片對應點對
4. camera calibration               相機成像參數
5. camera motion estimation         相機方位
6. surface point cloud generation   生成物體位置(點雲)
7. surface point cloud registration 對齊物體位置(點雲)
8. surface reconstruction           重建物體模型

Surface Point Cloud Generation

點雲生成。拍攝多張相片/一段影片,建立物體模型(點雲)。

triangulation:三角測量求得物體位置。使用兩張圖片甚至多張圖片,分別稱作two-view triangulation和multi-view triangulation。圖片座標是整數,兩條射線無法交於一點,因而衍生多種演算法。

depth estimation:物體深度是相機位置到物體位置的距離。藉由射線方向與物體深度可以輕易得到物體位置。使用兩張圖片甚至一張圖片,分別稱作binocular depth estimation和monocular depth estimation。大家從圖片抽取各種資訊,輔助判斷深度,例如邊緣偵測,因而衍生多種演算法。

bundle adjustment:同時求得相機方位、物體位置,一口氣完成兩個步驟。目前效果最好的方法。我認為原因是避開了邪教。

Triangulation

三角測量。一組對應點對,求得物體位置。

midpoint triangulation:公式解。效果最差。兩條射線的公垂線的中點,作為答案。

nonlinear triangulation:nonlinear least squares。效果最好。image feature matching的觀測值,camera imaging的估計值,兩點距離平方盡量小。此誤差稱作reprojection error

two-view triangulation:

min ‖y₁ - F(x; R₁,t₁,λ₁,f₁,⋯)‖² + ‖y₂ - F(x; R₂,t₂,λ₂,f₂,⋯)‖²
 x  

multi-view triangulation:
     N
min sum ‖yₙ - F(x; Rₙ,tₙ,λₙ,fₙ,⋯)‖²
 x  n=1

Hartley–Sturm linear triangulation:homogeneous linear least squares。效果尚可。一、camera imaging的估計值,採用y = λ K [R|t] x,省略λ = 1/zꟲ。導致估計值不在圖片上面。導致估計值從點座標變成射線向量。二、使用共線製造等式。令觀測向量和估計向量盡量共線(叉積盡量為零)。三、以焦距f求得點座標。

Depth Estimation

深度估計。一個像素,求得相機與物體的距離。

立體相機可以大幅簡化演算法,請見後面章節stereo camera之disparity estimation。光達可以直接得到深度,唉,轉來睏矣。

dense tracking and mapping (DTAM):multiobjective optimization。模仿disparity estimation的演算法。對應epipolar line,一口氣找到所有對應像素。

plane sweeping:image alignment。一、圖片座標:以第一張圖片為準,依序求得每個像素的深度。針對一個像素,枚舉深度,從中找到正確深度。二、對應座標:針對一種深度,推導物體座標,推導其餘N-1個圖片座標(三角測量)、像素數值(線性內插)。三、掃描面:以第一張圖片為準,當多個像素同時枚舉深度,則多個物體座標同時位於掃描面。四、圖片區域:針對一個像素,取得周遭圖片座標,形成圖片區域(自訂造型)。五、對應區域:針對一種深度,推導物體區域,推導其餘N-1個圖片區域(射影變換)、像素數值(線性內插)。六、對應區域差異:本應累計N(N-1)/2對差異,改為追蹤N-1對差異。每張圖片只跟第一張圖片比較差異。七、兩個區域差異:三種公式sum of absolute differences、sum of squared differences、zero-mean normalized cross correlation。採用後者。八、深度範圍:N-1對差異均低於臨界值者,推定為深度。否則視作未定義數值。九、掃描面方向:枚舉掃描面方向,迎合物體表面方向,減少未定義數值。

Bundle Adjustment

光束調整、光束平差。同時求得相機方位、物體位置。

bundle是指「束成一捆」。bundle adjustment是指「多條射線一併調整」。

calibration和adjustment互相呼應。calibration是指「讓儀器測量值可信可靠」,adjustment是指「讓儀器測量值迎合現況」。

大家採用最佳化。N個相機方位、M個物體位置。M組對應點群,一組有N個點。盡量減少reprojection error

          N   M
   min   sum sum ‖yₙᵢ - F(xᵢ; Rₙ,tₙ,λₙ,fₙ,⋯)‖²
Rₙ,tₙ,xᵢ n=1 i=1

做個總結:

上策同時求得相機方位與物體位置。下策先求相機方位,再求物體位置。上策採用最佳化演算法,下策用來估計疊代法初始值。

下策細分兩種方式。主角三角測量,直接求得物體位置。配角深度估計,先求物體深度,間接求得物體位置。主角用來計算答案,配角用來檢查答案。

此策略乃是針對相機。相機換成光達,策略地位上下翻轉。

Surface Point Cloud Registration

點雲對位。拍攝多張相片/一段影片,建立各種視角的物體模型(點雲)。拼裝物體模型(點雲),整合到世界座標系。

多張相片:大家採用iterative closest point

一段影片:大家自訂規則,調整數量與位置。【待補文字】

Surface Reconstruction

重建。拍攝多張相片/一段影片,重建物體模型。

點雲改成網格/參數曲面/貝茲曲面/體素/等值曲面/距離場。

網格/參數曲面/貝茲曲面:我沒有研究。【待補文字】

體素/等值曲面/距離場:早期作法是Poisson surface reconstruction。近期作法是deep learning,例如NeuralRecon

Scene Segmentation

分割。區隔每個物體。

早期作法是圖片做image segmentation,切割不同區塊。近期作法是每張圖片求得邊緣、深度、……,然後做deep learning。

Scene Graph Construction

建圖。聯繫每個物體,建立碰觸關係圖。

早期作法是兩兩圖片做image feature matching,連繫相似圖片。近期作法是每張圖片求得邊緣、深度、相機方位、……,然後做deep learning。

https://github.com/MIT-SPARK/Kimera
https://github.com/StanfordVL/3DSceneGraph

Scene Compositing(Match Move)

合成。拍攝多張相片/一段影片,將多餘物品塞入場景,配合相機方位。

Scene Synthesis(Scene Generation)

生成。家居擺設。

Motion Tracking🚧

簡介

Object Motion Estimation
(Object Motion Reconstruction)

運動估計。找到物體速度(位移速度、旋轉速度)。

我還不知道要怎麼做。

目前最精準的測量裝置是加速規陀螺儀,兩者合併為慣性測量單元實際範例實際範例

ego-motion = camera motion estimation in static scene
Scene Flow Estimation = optical flow with depth map

Moving Object Detection by 3D Flow Field Analysis
https://hal.science/hal-03565160/documentE
BundleSDF
https://bundlesdf.github.io/

Solving these problems will unlock
a wide range of applications in areas such as
augmented reality [34], robotic manipulation [22, 70],
learning-from-demonstration [71], and sim-to-real transfer.

Object Trajectory Estimation
(Object Trajectory Reconstruction)

軌跡估計。找到一連串物體位置(位移軌跡、旋轉軌跡)。

我還不知道要怎麼做。

[靜態物體]一、找出靜態物體(手工標記多個點)
      二、追蹤其圖片座標(image feature tracking)
      三、求得相機參數(camera calibration)
      四、求得靜態物體位置、相機位置(bundle adjustment)
[動態物體]五、找出動態物體(手工標記多個點)
      六、追蹤其圖片座標(image feature tracking)
      七、計算動態物體軌跡。

Object Reconstruction(Structure From Motion)

重建。deep learning硬幹?

Dynamic SfM / Dynamic SLAM
有時候是指"即時"偵測靜態場景,而不是偵測動態物體。

https://github.com/zhuhu00/Awesome_Dynamic_SLAM

Object Compositing(Match Move)

合成。人工調整模型位置?

Texture Tracking🚧

簡介

大家正在研究當中。章節架構無法定型。內容僅供參考。

Intrinsic Image

內在圖片。真相。去除光線影響,還原物體的真實顏色。

Image-based Modeling

拍攝相片,建立模型。可以想成是3D graphics的反運算。

平面變立體、2D變3D。根據像素顏色,推論模型形狀。像素值,是素材與光線的綜合結果。去除素材影響,則得光線。再從光線強度,推估光線入射角,得到表面法向量。

這件事聽起來很蠢。想要建立模型,直接使用偵測模型的器材即可,何必使用相機?原因是幾乎人人隨身攜帶智慧型手機,那上面裝著計算機和相機;但是不是人人都有深度感測儀。

Image-based Texturing

貼材質。在模型上貼材質。

Image-based Relighting(Image-based Rendering)

重新打光。偵測光源與光強,然後重新打光。

枚舉光源位置:架設大量光源,紀錄物體呈現的顏色。以內插得到特定光源位置的顏色。架設大量攝影機,以便重建模型。

Light Tracking🚧

簡介

大家正在研究當中。章節架構無法定型。內容僅供參考。

光線的數學知識:Radiance

radiance irradiance luminance illuminance。

https://computergraphics.stackexchange.com/questions/7503/

光源的亮度、物體表面特定方向的亮度。

Environment Light Estimation

環境光源估計。偵測來自各種方向360度的光強。

拍攝金屬球:相機對準金屬球中心,拍攝金屬球的半面。根據相片的金屬球的半面的表面亮度,運用光線入射角等於折射角的原理,推測半個環境的光源與光強。調整相機方位(例如水平周轉一圈),以各種方位拍攝,得到整個環境的光源與光強。

拍攝平面鏡:相機拍攝平坦的反光物體,以各種方位拍攝(例如直線平移),推測光源與光強。

http://www.johanneskopf.de/publications/gdibr/
https://snsinha.github.io/proj/ReflectionsIBR/

Light Source Estimation / Lighting Estimation

光源估計。偵測光源位置暨方向。

打光估計。直接輸出合成場景打光之後的照片。

光源可能是點光源、有向光源、面光源、環境光源、……。

早期作法是運用Graphics領域的光線追蹤演算法與照明模型,推測光源與光強。稱作光源估計。

https://www3.cs.stonybrook.edu/~cvl/illuminationModeling.html

近期作法是事先找到圖片高光區塊,取樣,Gaussian,deep learning催落去萬事OK。稱作打光估計。

ElasticFusion
https://github.com/mp3guy/ElasticFusion
DeepLight
https://www.cg.tuwien.ac.at/research/publications/2019/kan-2019-dli/
EMLight / GMLight
https://github.com/fnzhan/EMLight

Radiance Field Estimation

照明場估計。找到空間每一處的光線亮度。

neural radiance field (NeRF):deep learning。光場變體素。再以voxel rendering繪製物體模型。光場相機的發明人。

Gaussian splatting:deep learning。彩色點雲變圖片。彩色點雲,套用高斯分布(外觀宛如水花噴濺),以便平滑合成。位置平滑合成是blobby surface。顏色平滑合成是subsurface scattering。

Photo Stitching🚧

簡介

Image Straightening

矯直。調整相機方位,擺正相片裡面的物體。

簡單來說,相機做剛體變換(圖片內容做射影變換),使得相機座標軸對齊物體座標軸(圖片座標軸對齊物體座標軸)。

進階應用有變形圖片GPS

矯直方式通常是讓使用者手動劃記直線所在位置,然後自動實施image warping。

Image Rectification

矯正。調整相機方位,使得相機正面方向一致、相機正面垂直於相機連線(圖片共面、圖片橫軸平行於相機連線)。

簡單來說,相機做剛體變換(圖片內容做射影變換),使得相機座標軸相互對齊(圖片座標軸相互對齊)。

事先實施essential matrix estimation求得相機相對方位。

矯正方式有無限多種。追加約束條件,以找到你喜歡的那一種。

固定相機壹方位:相機貳從該處橫向位移,各種長度都是答案。

固定相對位移:相機連線作為轉軸,各種角度都是答案。

又矯直又矯正:我沒有研究。

Image Stitching(Image Mosaicing)(Panorama)

縫合。拍攝多張相片,源自各種方位。銜接成一張寬廣相片,稱作環景照。經典範例是谷歌街景教學文章程式碼

每張相片因為光線角度導致色彩差異、因為廣角鏡頭導致景物扭曲。需要校正亮度、扳直曲線、填補像素、裁切圖片。

Image Stabilization🚧

簡介

Video Stabilization

安定。俗稱防手震。拍攝影片,相機震動搖晃,相機位置不斷改變,導致像素位移。令影片當中每張圖片內容互相對齊。

目前多以物理方式處置。花14美金自己DIY

現成影片只好以數學方式處置。optical flow之類的

Video Magnification

放大。以高幀率高解析度的攝影機進行拍攝,觀察物體極短時間、極小範圍的微小變化。應用範圍很廣

Stereo Camera🚧

簡介

相機的數學知識:Projective Geometry

射影幾何。一臺相機,拍攝一個物體,得到一張相片。

只從圖片無法得知物體遠近。一個圖片座標對應無限多個物體座標。

相機的數學知識:Epipolar Geometry

對極幾何。兩臺相機,拍攝一個物體,得到兩張相片。主題曲

epi-是指「在此之上、額外」,我不知道為何翻譯成「對」。

兩張圖片可以判別物體遠近。一對圖片座標對應一個物體座標。

對極幾何的元件名稱。

baseline        基準線 兩臺相機所在位置連線
epipole         對極  圖片平面與相機連線的交點
epipolar plane  對極面 三角測量平面
epipolar line   對極線 圖片平面與三角測量平面的交線

對極幾何的元件性質。

一、沿著物體表面移動。
  相機方位、基準線、圖片方位、對極,固定不變。
  物體座標、對極面、圖片座標、對極線,不斷改變。
二、一張圖片的各種對極線,通通穿過同一個對極。
三、相機一前一後,保證讓對極出現在圖片上。
一、兩張圖片的對極的景物一致(如果景物沒被遮擋)。
二、兩張圖片的對應點的景物一致(如果景物沒被遮擋)。
三、兩張圖片的對極線的景物一致(如果景物沒被遮擋)。
  沿著物體表面移動。
  令物體座標逐漸改變、不脫離對極面;
  則圖片座標隨之改變、不脫離對極線。
  注意到,間距疏密不同。

Stereo Camera(Stereoscopy)

立體相機。兩臺相同相機,已經矯正方位。

相機矯正,對極幾何的元件亦隨之矯正。

相機正面垂直於baseline。
圖片平行於baseline。
epipole在無窮遠處。
一張圖片的epipolar line都是水平線。
兩張圖片的對應epipolar line是相同水平線。

epipolar line = scanline。得以精簡計算過程。

Disparity Estimation

視差估計。找到每個像素的視差。

對應點對,依照其圖片座標,描繪於同一張圖片,左圖或右圖皆可。視差是此點對的水平偏移距離(無窮遠物體,視差為零)。

視差也是對應點對到圖片中央的水平距離加總(當心正負號)。

baseline length : depth = disparity : focal length。已知相機方位和焦距,藉由視差可得深度。進一步得到點雲。

block matching:dynamic programming。一、最佳匹配:對應scanline,一口氣找到所有對應像素。對應像素位置單調遞增。左右視野看見不同部件必須刪除。綜上所述就是longest common subsequencedynamic time warping。二、匹配成本:左圖每個像素的每種視差的成本C(p,d)。左圖像素p利用視差d取得右圖對應像素p',考慮兩個像素的差異。三、像素差異:針對一個像素,取得周遭像素,形成方形區塊(自訂大小),考慮兩個區塊的差異。四、區塊差異:三種公式sum of absolute differences、sum of squared differences、zero-mean normalized cross correlation。採用後者。五、狀態轉移成本:視差差異(乘上自訂權重)。六、匹配範圍:事先做image segmentation。景色陡變之處通常是視差陡變之處。

plane sweeping:同上,對調內外迴圈,先枚舉d再枚舉p,宛如掃描面。就這樣。

semi-global matching:dynamic programming。一、成本:左圖每個像素的每種視差的成本C(p,d)。左圖像素p利用視差d取得右圖對應像素ebm(p,d)。對應像素的差異,作為成本。原論文採用複雜的cross-entropy,MATLAB採用簡潔的census transform的Hamming distance。二、追加成本:視差差一的鄰居像素成本P1、視差差二以上的鄰居像素成本P2。令P1 < P2。令P2是P2'除以對應像素數值差距,差距越多P2越小。三、修改成本公式:鄰居改成16方位。一般公式改成遞迴公式(方形區塊改成任意造型)。一個方位的成本Lr(p,d)。所有方位的成本S(p,d)。四、視差:針對一個像素p,令成本S(p,d)最小的視差d,作為答案。五、cross check test:再以右圖為主角,計算答案。答案差異不可超過一,否則視作未定義數值。

View Synthesis

視野生成。根據特定相機方位的相片,產生其他相機方位的相片。

structure from motion是循規蹈矩的作法,需要大量相機方位,以便建立模型、重新成像。view synthesis是投機取巧的作法,只要兩三個相機方位,不去建立模型,而是改造相片。

Scharstein's algorithm:矯正兩張圖片。各自求得視差。各自求得新圖片。對齊兩張新圖片。合成與補洞。反向矯正新圖片。

Light Field Camera🚧

簡介

相機的數學知識:Light Field

光場。三維空間的每個地點、每種方向,各有一個光線顏色數值。空氣所在之處是透明、不存在數值。

Light Field Camera

光場相機。以相機陣列篩選光線方向,以照片陣列推導光場。

camera array:大量相機,排列整齊,形成二維陣列平面。藉由鏡頭位置,篩選光線方向。例如台灣特產方格壓花玻璃窗

plenoptic camera:合併鏡片,保留鏡孔。換句話說,大量針孔相機,排列整齊,形成二維陣列平面,外頭共用一個鏡頭。

Light Field Display

光場顯示器。光線路徑反過來。可以顯示立體影像。

Light Field Manipulation

光場操作。利用光場,模擬各種光學機制、成像機制。

例如重新對焦,得到各種焦距的相片。例如更換鏡頭,得到各種視野的相片。然而普通相機也能達成相同效果。等你想出殺手級應用。