Photometric Stereo 光度立体三维重建(四)——光源标定
光源標定是進行光度立體三維重建的第一步,本文將介紹兩種光源標定方法——基于金屬球反射的標定以及基于“SFM”思想的標定
1.基于金屬球反射的標定
標定光源的一種方法是使用金屬球,在排到的金屬球的照片上面的最亮的點指明了光源的方向
來源于http://pages.cs.wisc.edu/~csverma/CS766_09/Stereo/stereo.html的示意圖:
然而,這幅示意圖的幾何向量標注具有一定的誤導性,我將其換個角度表示以便展開光源方向的推導:
其中,反射方向R是固定的(因為它一直正對著觀察者):
R=[0,0,1]T\ \mathbf{R}=[0,0,1]^T?R=[0,0,1]T
我們取Px,Py\ P_x,P_y?Px?,Py?為最亮點在圖像中的像素坐標位置,Cx,Cy\ C_x, C_y?Cx?,Cy?為金屬球在圖像中的中心點(可以由mask圖像確定),由此我們可以求出法線向量N=[Nx,Ny.Nz]T\ \mathbf{N}=[N_x,N_y.N_z]^T?N=[Nx?,Ny?.Nz?]T:
Nx=Px?Cx\ N_x=P_x-C_x?Nx?=Px??Cx? Ny=Py?Cy\ N_y=P_y-C_y?Ny?=Py??Cy? Nz=(R2?Nx2?Ny2)\ N_z = \sqrt{(R^2-N_x^2-N_y^2)}?Nz?=(R2?Nx2??Ny2?)?
N_z中的R是球的半徑
求得法線向量后,我們可以知道:
R+L=N=2S\ \mathbf{R}+\mathbf{L}=\mathbf{N}=2\mathbf{S}?R+L=N=2S L=2S?R\ \mathbf{L}=2\mathbf{S}-\mathbf{R}?L=2S?R
我們現在轉移到求S,因為S是R在N上的投影(向量投影計算可參考這里),因此可以得到:
S=(R?N)∣N∣2N=(R?N)N\ \mathbf{S}=\frac{(\mathbf{R}·\mathbf{N})}{|\mathbf{N}|^2}\mathbf{N}=(\mathbf{R}·\mathbf{N})\mathbf{N}?S=∣N∣2(R?N)?N=(R?N)N
帶入L的表達式,可得:
L=2(N?R)R?R\ \mathbf{L}=2(\mathbf{N·\mathbf{R}})\mathbf{R}-\mathbf{R}?L=2(N?R)R?R
我們便可以利用上式求解光源方向
2.基于“SFM”的光源標定
論文標題:Light Structure from Pin Motion: Geometric Point Light Source Calibration(ECCV2018 / IJCV2020)
2.0 SFM簡單回顧
SFM(structure from motion)目的是從圖像特征點的對應同時恢復3D結構(3D點的位置)和姿態(相機的外參數)
設Pw\ P_w?Pw?是世界坐標系的一個點,Pc\ P_c?Pc?是Pw\ P_w?Pw?在相機坐標系中的位置,則有Pc=[Rt]Pw\ P_c=[R \quad t]P_w?Pc?=[Rt]Pw?,該點在成像平面的齊次坐標為p=[μ,v,1]T\ p=[\mu,v,1]^T?p=[μ,v,1]T,其在歸一化平面的坐標為x=[Xc/Zc,Yc/Zc,1]\ x=[X_c/Z_c,Y_c/Z_c,1]?x=[Xc?/Zc?,Yc?/Zc?,1],有p=Kx,x=K?1p\ p=Kx,x=K^{-1}p?p=Kx,x=K?1p,則通過對極幾何約束,我們可以獲得同一個特征點在兩幅圖像(相機在兩個姿態)下坐標的對應關系:
對于本質矩陣
E=[t]×R\ E=[t]_\times R?E=[t]×?R
和基礎矩陣
F=K?TEK?1\ F=K^{-T}EK^{-1}?F=K?TEK?1
本質矩陣表征了點在兩個姿態下歸一化平面坐標的對應關系:
x2TEx1=0\ x_2^TEx_1=0?x2T?Ex1?=0
基礎矩陣表征了點在兩個姿態下齊次像素坐標的對應關系:
p2TFp1=0\ p_2^TFp_1=0?p2T?Fp1?=0
考慮E\ E?E的尺度等價性,我們可以使用八點法來求解基礎矩陣E\ E?E
在求出E\ E?E后,可以通過SVD分解求出R和t\ R和t?R和t
如果場景中的特征點都落在同一平面上,則可以通過單應矩陣H\ H?H來進行運動估計,單應矩陣描述了兩個平面間的映射關系,對于圖像I1,I2\ I_1,I_2?I1?,I2?中匹配好的落在同一個平面上的特征點p1,p2\ p_1,p_2?p1?,p2?,其對應關系為:
p2=K(R?tnTd)K?1p1=Hp1\ p_2=K(R-\frac{tn^T}ze8trgl8bvbq)K^{-1}p_1=Hp_1?p2?=K(R?dtnT?)K?1p1?=Hp1?
單應矩陣可以通過4對匹配特征點算出(特征點不能三點共線)
在估計出相機運動后,可以通過三角測量估計出地圖點的深度。
2.1 方法介紹
在第一節所介紹的用金屬球進行標定的方法中,存在以下的弊端:
1.需要擁有相對精確的球體
2.需要注釋出每個圖像的高亮中心
3.需要標出圖像中球的輪廓(mask)
4.小的標注誤差會導致光源位置大的偏差
因此作者提出一種新的方法來進行光源的標定:作者通過在標定板上插上帶有一定大小頭部的針,通過在固定相機下改變標定板的姿態,從而獲得當前光源的位置信息:
論文提出方法的流程圖:
2.1.1陰影的形成模型
假設陰影接收平面π固定于世界坐標系的x-y平面
一個近點光源在世界坐標系的位置為:
l=[lx,ly,lz]T∈R3\ \mathbf{l}=[l_x,l_y,l_z]^T \in \mathbb{R}^3?l=[lx?,ly?,lz?]T∈R3
假設陰影制造物(caster),即針頭的位置為c∈R3\ \mathbf{c} \in \mathbb{R}^3?c∈R3,它產生的陰影在π平面上位置為s∈R2\ \mathbf{s} \in \mathbb{R}^2?s∈R2,在世界坐標系中位置是s ̄=[sT,0]T\ \overline{\mathbf{s}}=[\mathbf{s}^T,0]^T?s=[sT,0]T,因為l,c,s ̄\ \mathbf{l},\mathbf{c},\overline{\mathbf{s}}?l,c,s在同一直線上,因此由其中兩個點組成的任意兩直線是平行的,因此可以得出:
將
代入(1)式得:
將叉積展開得到:
從中可以求出sx和sy:
我們將s寫為齊次坐標,并使用縮放因子γ\ \gamma?γ,得到:
繼續整理得:
將等號右邊改寫為矩陣相乘的形式得:
我們稱L\ \mathbf{L}?L為光矩陣,我們將其分解可以得到其內參和外參:
以上是關于近點光源的,對于遠點光源,光線是平行的,因此光向量應該指光源的方向,而不是光源點的位置,因此(1)變為:
同理我們可以整理為:
我們發現對于遠點光源,光矩陣(3,3)的位置是0,這個光矩陣類似于在相機矩陣中的正交投影
我們希望將近點光源和遠點光源的光矩陣統一起來,在齊次坐標中,相差一個常數的縮放是等價的,因此我們把兩個光矩陣都除以lz\ l_z?lz?,得到:
可以發現,當光源無限遠時,近點光源光矩陣變為遠點光源光矩陣,因此我們使用下式作為統一的陰影投影等式:
這也是為什么論文名為Light Structure from Pin Motion的原因,論文中恢復光源方向的思想是跟SFM類似的,點光源和針孔相機可以被相似的數學模型所描述,其相似的對應關系為:
1.點光源=針孔相機
2.陰影接收面=圖像平面(成像面)
3.陰影caster=場景點
4.光矩陣L=相機投影矩陣P=KT=K[R∣t]\ P=KT=K[R|t]?P=KT=K[R∣t]
Light Structure from Pin Motion與SFM的對比:
2.1.2 陰影對應中的對極幾何
回顧SFM,相機運動的求解是基于對極幾何約束的,那么,在Light Structure from Pin Motion中,也應該找出陰影對應的對極幾何約束。
作者移動光源,同時casters和陰影接收面是靜止的,在其中分析陰影的對應關系,這類似于SFM中移動的相機與靜置的場景點。
陰影的對極幾何關系圖:
caster是靜置且位置未知的,對于由光源l1\ l_1?l1?投影得到的陰影s\ s?s,caster可以在線段l1s\ l_1s?l1?s上的任意位置,當光源位置移動到l2\ l_2?l2?時,一系列陰影的可能位置形成一條線,這相當于相機對極幾何中的極線。
我們現在希望構造出陰影對極幾何中的基礎矩陣F\ F?F
記c\ c?c是caster的位置,l1≠l2\ l_1 \neq l_2?l1??=l2?是在標定目標坐標系下的光的位置,Li是li\ L_i是l_i?Li?是li?的光矩陣,L1+=(L1TL1)?1L1T是L1\ L_1^+=(L_1^TL_1)^{-1}L_1^T是L_1?L1+?=(L1T?L1?)?1L1T?是L1?的偽逆,?L1∈null(L1)\ \empty_{L_1} \in null(L_1)??L1??∈null(L1?)是L1\ L_1?L1?一維零空間下的非零向量,η\ \eta?η是一個縮放常量,我們可以得到:
對上式左乘s~2T[L2?L1]×\ \tilde{s}_2^T[L_2 \empty_{L_1}]_\times?s~2T?[L2??L1??]×?,得到:
由此我們得到了表征陰影對應的基礎矩陣,對于Li\ L_i?Li?
我們得到零空間向量?L1\ \empty_{L_1}??L1??
將其代入(5)式得到基礎矩陣:
這個矩陣是反對稱的,對于陰影s1~=[u,v,1]T,s2=[u′,v′,1]~\ \tilde{s_1}=[u,v,1]^T, \tilde{{s_2}=[u',v',1]}?s1?~?=[u,v,1]T,s2?=[u′,v′,1]~?,陰影的對極約束為:
因此我們可以解齊次線性方程組來求解在一個縮放因子下成立的系數f:
這實際上相當于估計只進行純平移而不進行相對旋轉的相機的常規的本質矩陣,由式(2)可以看出,旋轉矩陣為單位陣
2.2 模型求解
論文提出的方法中,通過在固定點光源下從固定視點(相機位置固定)對標定目標進行多次觀察,同時改變標定目標的姿態,自動實現點光源標定。陰影caster相對于校準目標的3D位置是未知的,這使得它特別容易建立目標,而問題仍然是易于處理的。
2.2.1 通過光束平差法標定光源方向
目標:通過觀察投影自未知caster的陰影求出(4)中的l\ l?l。一次觀察不足以提供充分的信息,因為我們讓陰影接收面經歷多個姿態{[Ri∣ti]}\ \{[R_i|t_i]\}?{[Ri?∣ti?]}
在第i個姿態中,光源在接收面坐標系中的位置li\ l_i?li?與在世界坐標系中的位置l\ l?l的關系為:
光矩陣Li\ L_i?Li?為:
作者不僅使用多個姿態,更使用多個casters{cj}\ \{c_j\}?{cj?}來提高校準精度,因此對于姿態i和caster j,我們得到陰影sij\ s_{ij}?sij?,此時(4)變為:
假設目標姿態{[Ri∣ti]}\ \{[R_i|t_i]\}?{[Ri?∣ti?]}已知,我們目標是估計光源在世界坐標系中的位置l\ l?l以及caster在標定目標坐標系中的位置cj\ c_j?cj?,我們通過最小化重投影誤差構建一個最小二乘問題:
可以使用LM算法來求解這個非線性最小二乘問題,為了魯棒的估計使用隨機抽樣一致性(RANSAC)算法:我們不斷選取一個隨機的觀測集,估計(l,cj,λij)\ (l,c_j,\lambda_{ij})?(l,cj?,λij?),選擇有最小殘差的一個
2.2.2 Bundle Adjustment的初始化
(1)近點光源
類似于式(1),我們可以寫出:
將
代入得到:
其中,
我們可以重寫上式為:
擴展叉積得到:
我們將它整理一下,得到:
其矩陣形式為:
這個方程包含一個觀測結果,即姿態i和caster j的結合,但是我們需要堆砌多個觀測下的方程。為了簡化步驟,我們先將(18)分為幾個子矩陣:
讓Np,Nc\ N_p,N_c?Np?,Nc?為目標姿態和casters的數量,可以得到整個系統的方程為:
因為所有的觀測共享l=[lx,ly,lz]T\ l=[l_x,l_y,l_z]^T?l=[lx?,ly?,lz?]T,每個θj\ \theta_j?θj?有12個未知數,因此我們有3+12Nc\ 3+12N_c?3+12Nc?個未知數
為了防止外點的影響,我們通過L1最小化來魯棒求解:
為了求解此式,我們要滿足:
因此,不管casters有多少,5個姿態下的觀測足以求解方程
解出θ?\ \theta^*?θ?后,我們忽視二階變量如lxcj,x\ l_xc_{j,x}?lx?cj,x?,使用cj?,l?\ c_j^*,l^*?cj??,l?來初始化(17)
(2)遠點光源
對于遠點光源,矩陣A的秩為3+12Nc?1\ 3+12N_c-1?3+12Nc??1,因此,不幸的是,我們不能使用(18)來處理遠光,但我們可以自動檢測這種情況,切換到對遠光的方程,求解,并切換回(17)的BA。因此,用戶不必為他們的場景選擇一個光照模型。
我們通過構造近點光源的矩陣A,然后檢測A的奇異值,如果最大和最小的奇異值之比大于4e+4,我們轉換到遠光求解
對于遠光,類似于(3),我們使用li=RiTl\ l_i=R_i^Tl?li?=RiT?l重寫為:
類似得,我們將其整理得:
設置l=[lx,ly,1]T\ l=[l_x,l_y,1]^T?l=[lx?,ly?,1]T(2自由度)得到:
整理得:
將其改寫為矩陣形式,得:
其子矩陣的形式為:
多次觀測的方程堆疊和求解類似于(19)(20),為了求解這個線性方程組,我們要滿足:
因此只需要4個姿態下的觀測,因為對于遠光的方程組求解得到的是光的方向(direction),但光束平差法需要光的位置(position),因此我們需要它們的一個轉換。我們從其中一個casters開始,將光源不斷移動到很遠的地方:
c是世界坐標系下任意的caster的位置,k是大常數(10e+10),hc是caster的高度
2.2.3 陰影對應
在(17)(18)(21)中,我們需要得到在不同圖像Ii\ \mathbf{I}_i?Ii?中同一個caster j產生陰影sij\ s_{ij}?sij?的對應。形式上,陰影對應搜索意味著我們需要找到與圖像之間對應陰影匹配的置換
設S\ S?S為陰影{si~}\ \{\tilde{s_i}\}?{si?~?},S′\ S'?S′為陰影{si′~}\ \{\tilde{s'_i}\}?{si′?~?},將其水平疊加成矩陣,對于基礎矩陣,我們尋找:
也就是說,S′\ S'?S′的每一列都是一個陰影向量,我們將置換矩陣P′\ P'?P′右乘,以對各個陰影向量進行排列,只有S′和S\ S'和S?S′和S的每一列對應的陰影向量是在兩幅圖像中一一對應時,基礎矩陣的對極約束等式才為0
不同于傳統的圖像特征匹配,我們不能使用特征描述符來縮小搜索范圍,因為我們希望陰影非常小(為了讓陰影的中心可以精確定位),所以不能改變caster的形狀以使其陰影與其他caster的陰影有明顯區別
作者通過檢測多幅圖像對應的一致性以獲得陰影對應
對應一致性檢驗
作者建立了兩個圖像池:已經建立了陰影對應的圖像——確定池(established pool) 以及對應未知的圖像——未知池(unknown pool)
確定池用一張隨機的,未知的圖像進行初始化,目的是將盡可能多的圖像從未知池移動到確定池中
工作分為兩個階段,在第一個階段中,在確定池隨機選取圖像Ie\ I_e?Ie?,在未知池選取k張圖像Ii1,...Iik\ I_{i_1},...I_{i_k}?Ii1??,...Iik??
假設
是一個二值函數,其為真的條件是:當且僅當對于圖像Ia,Ib\ I_a,I_b?Ia?,Ib?,Ia\ I_a?Ia?中的陰影si\ s_i?si?和Ib\ I_b?Ib?中的陰影sj\ s_j?sj?獲得匹配,使得式(22)最小化
然后,如果:
成立(即,從確定池抽取的圖像開始,與未知池的第一張圖像匹配,然后未知池的第一張與第二張圖像進行匹配,以此類推,通過第一張(確定池圖像)到最后一張(未知池的最后一張圖像)能沿著圖像鏈一直對應下去),我們就將未知池的圖像Ii1,...Iik\ I_{i_1},...I_{i_k}?Ii1??,...Iik??移動到確定池,為了使約束嚴格,作者在文章中選擇k=3,一旦有一半的圖像轉移到了確定池,我們就轉到第二個階段
在第二個階段中,假設確定池中所有的圖像都是一致對應的。因此,如果我們考慮一個未知圖像和一個caster,那么所有已確定的圖像中對應于那個特定caster的所有陰影都應該匹配未知圖像中的相同陰影,這適用于所有的陰影caster。我們隨機選取一幅未知圖像,驗證這一準則,如果超過一半的已確定圖像與未知圖像的對應關系一致,將其移動到確定池,否則將其丟棄,當未知池為空時,第二階段結束
陰影檢測
作者使用模板匹配方法進行陰影檢測
作者生成了由一條線和一個圓組成的陰影合成圖像作為模板。為了處理變化的射影變換,生成的模板具有12個旋轉角度,每個角度都有3個大小的模板(也就是有36個模板)。作者將輸入圖像二值化后再進行模板匹配。進一步,作者使用針頭頭部的顏色(紅色)來區分頭部和頭部陰影
2.3 實際操作(代碼)
論文的github倉庫:https://github.com/hiroaki-santo/light-structure-from-pin-motion
(1)相機標定
論文的代碼中是使用Arcuo碼進行標記的,因此我們首先要制作Arcuo標定板:
打印得到的Aruco標定板:
我在測試時使用的是手機攝像頭,通過不同角度拍攝屏幕上顯示的標定板,獲得了14張圖片來進行標定:
標定得到的數據:
上面是手機相機的內參矩陣,下面是5個失真系數
代碼里的detect_markers.py可以檢測aruco標記:
剩下的陰影檢測,實際光源標定的步驟在代碼的README文檔里也有操作指導,這里就不展開了。
我組建了一個光度立體技術的交流群,有興趣的朋友可以一起來討論一下!
總結
以上是生活随笔為你收集整理的Photometric Stereo 光度立体三维重建(四)——光源标定的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Java 同步器
- 下一篇: 矩阵位移法matlab编程,矩阵位移法_