あらいさん日記（仮）

あらい日記

活動の記録

Spinnaker SDKでFLIRのカメラをPythonから動かす

機材

FLIRのSpinnaker SDKではC/C++,Pythonといった主要なプログラミング言語をサポートしています．ここではPythonを使ってFLIRのカメラを動かすための具体的な方法を説明します．

Spinnaker SDKのインストール方法はこちらを参考にしてください．

PySpinのインストール

Spinnaker SDKをPythonから使うためには，PySpinをインストールする必要があります．
numpyとmatplotlibをインストールしてない場合は，pipでインストールします．Pillowは無くても動きますが，あった方が良いみたいです．

MacOS

  1. /Applications/Spinnaker/PySpin内にあるspinnaker_python-2.0.0.109-cp37-cp37m-macosx_10_9_x86_64.tar.gzを解凍します．2.0.0.109はSpinnaker SDKのバージョン，cp37はPythonのバージョンです．
  2. 解凍したフォルダ内にあるPySpinのwheelをインストールします．
    pip install spinnaker_python-2.0.0.109-cp37-cp37m-macosx_10_9_x86_64.whl

PySpinの使用例

FLIRが配布しているコードの中には，PySpinのexamplesがありますが，なぜか非常に複雑です．ここではPySpinを使うためのコードを要点を絞って紹介します．

撮影・保存

import PySpin

#カメラの指定・初期化
system = PySpin.System.GetInstance()
cam_list = system.GetCameras()
cam = cam_list.GetByIndex(0)
cam.Init()

#撮影開始
cam.BeginAcquisition()

#カメラから画像を転送＆保存
image_result = cam.GetNextImage()
image_converted = image_result.Convert(PySpin.PixelFormat_Mono8, PySpin.HQ_LINEAR)
image_converted.Save("image.png")
image_result.Release()

#撮影終了
cam.EndAcquisition()

#終了処理
del cam
cam_list.Clear()
system.ReleaseInstance()

注意したいのは終了処理（特にdel camやcam_list.Clear()）の部分で，これが実行されない場合エラーになります．

numpy配列への変換

image_NDArray = image_result.GetNDArray()

PySpinでは撮影した画像をnumpy配列に変換することができます．そのため，opencvやnumpyとの連携が容易です．

撮影の設定

カメラのパラメータを設定をするにはQuickSpinAPIかGenAPI(GenICamAPI)を使います．QuickSpinAPIは，ほとんどのパラメータがサポートされていて，GenAPIより短いコードで実行できます．そのため，基本的にはQuickSpinAPIを使えば良いです．

f:id:raccoon15:20200430093255p:plain
設定できるパラメータはSpinViewの左下のメニュー(Features)から確認できる
ここではQuickSpinAPIを使った各種パラメータの設定方法を紹介します．

露光時間

#Auto exposureをon (continuousモード)
cam.ExposureAuto.SetValue(PySpin.ExposureAuto_Continuous)

#Auto exposureをoff
cam.ExposureAuto.SetValue(PySpin.ExposureAuto_Off)
#Exposure modeをTimedに設定
cam.ExposureMode.SetValue(PySpin.ExposureMode_Timed)
#Exposure timeを20000usに設定
cam.ExposureTime.SetValue(20000)

#設定されているexposure timeを取得
value = cam.ExposureTime.GetValue()

ゲイン

#Auto gainをon (continuousモード)
cam.GainAuto.SetValue(PySpin.GainAuto_Continuous)

#Auto gainをoff
cam.GainAuto.SetValue(PySpin.GainAuto_Off)
#Gainを10.5dBに設定
cam.Gain.SetValue(10.5)

#設定されているGainを取得
value = cam.Gain.GetValue()

ガンマ

#Gammaを1.5に設定
cam.Gamma.SetValue(1.5)

#設定されているGammaを取得
value = cam.Gamma.GetValue()

黒レベル

#BlackLevelSelectorをAllに設定
cam.BlackLevelSelector.SetValue(PySpin.BlackLevelSelector_All)
#BlackLevelを1.5%に設定
cam.BlackLevel.SetValue(1.5)

#設定されているBlackLevelを取得
value = cam.BlackLevel.GetValue()

転送する画像の順番

#GenAPI
#転送する画像を最新のもの(NewestOnly)に設定
s_node_map = self.cam.GetTLStreamNodeMap()
handling_mode = PySpin.CEnumerationPtr(s_node_map.GetNode('StreamBufferHandlingMode'))
handling_mode_entry = handling_mode.GetEntryByName('NewestOnly')
handling_mode.SetIntValue(handling_mode_entry.GetValue())

webカメラのように撮影した画像をリアルタイムで確認したい場合は，StreamBufferHandlingModeをNewestOnlyに設定して，最新の画像のみを転送してくるようにします．OldestFirstに設定すると，カメラ内のバッファに保存されている古い画像から順に転送されます．

OpenCVのVideoCaptureクラスっぽく使う

PySpinはカメラの詳細な設定をできますが，それ故に複雑で独特なコーディングになりがちです．一方で，OpenCVのVideoCaptureクラスは簡素なコードでwebカメラを動かすことができます．そこで，OpenCVのVideoCaptureクラスっぽいコードで，FLIRのカメラを動かすためのラッパーを作ってみました．複雑なコードでお困りの方は使ってみてください． github.com

FLIRのSpinnaker SDKのインストール方法のメモ

機材

FLIRのマシンビジョンカメラ（BlackFly, Oryx等）はSpinnaker SDKを使って動かすことができます．ここでは，Spinnaker　SDKをインストールして，FLIRのカメラを動かすための方法をメモしておきます．

Spinnakerのインストール

下のリンクからSpinnakerSDKをダウンロードします．すべてのバージョンをダウンロードしようとすると時間がかかるので，自分の環境に適したものを選択してダウンロードします． www.flir.com

MacOS

  1. Homebrewでpkg-config,libomp,libusbをインストールします．
    brew install pkg-config libomp libusb
  2. MacOS/Spinnaker-2.0.0.109.dmgを開き，Spinnaker-2.0.0.109.pkgから手順にしたがってインストールを進めます．
  3. /Applications/Spinnaker/apps/SpinView_QT.appを実行すると，カメラからの映像を確認することができます．

Windows

  1. SpinnakerSDK_WEB_1.29.0.5_x64.exeもしくはSpinnakerSDK_FULL_1.29.0.5_x64.exeを実行して，手順にしたがってインストールを進めます．
    Camera EvaluationかApplication Developmentを選択します．どちらを選んでもSpinView（GUIのビューワー）はインストールされ，後者のApplication DevelopmentはSDKが付属します．
  2. SpinViewを実行すると，カメラからの映像を確認することができます．

SpinViweの比較

SpinViewはFLIRのカメラの設定・撮影・保存といった動作をGUIで実行できます．

f:id:raccoon15:20200430130600p:plainf:id:raccoon15:20200430130604p:plain
SpinViewの比較（左）MacOS，（右）Windows
MacOS版とWindows版のSpinViewのスクリーンショットを比較すると，Windows版だけが持つ機能を確認できます．Windows版だけの機能は，「ヒストグラムの表示」「スライダーなどのUIによるパラメータの設定」「ログビュアー」などがあります．

偏光カメラで出来ること

偏光

概要

偏光を見ることで，通常のカメラでは得られない情報を獲得することができます．最近では偏光カメラが比較的安価に手に入るようになり，偏光イメージングの機運が高まっています．ここでは，偏光カメラで撮影した画像を使って，偏光イメージングの応用例を紹介します．

偏光の基礎

＊ここでは直線偏光のみを考えています．

偏光の状態 

偏光と偏光板

偏光板を回転させたときに観測される輝度Iは次式のようなコサイン関数で表せられます．

I(\upsilon) = \frac{I_{max}+I_{min}}{2} + \frac{I_{max}-I_{min}}{2} \cos(2\upsilon - 2\psi)

ここで\upsilonは偏光板の偏光角，I_{max}とI_{min}は輝度の最大値と最小値，\psiは位相角です．

f:id:raccoon15:20200103182658p:plain

偏光板を通して観測される光の輝度

偏光の度合いを表す偏光度\rhoは次式で表せられます．完全偏光の場合は\rho=1，ランダム偏光の場合は\rho=0になります．

\rho=\frac{I_{max}-I_{min}}{I_{max}+I_{min}}

ストークスベクトル

偏光の状態はストークスベクトルsを用いて表すことができます．

s = \begin{pmatrix}s_{0} \\ s_{1} \\ s_{2}\end{pmatrix}=\begin{pmatrix}I_{max}+I_{min} \\ (I_{max}-I_{min})\cos2\psi \\ (I_{max}-I_{min})\sin2\psi \end{pmatrix}=(I_{max}+I_{min})\begin{pmatrix}1 \\ \rho\cos2\psi \\ \rho\sin2\psi \end{pmatrix}

s_{0}：光の輝度の合計．
s_{1}：0^{\circ}方向の直線偏光成分．
s_{2}：+45^{\circ}方向の直線偏光成分．

コサイン関数I(\upsilon)のパラメータI_{max}，I_{min}，\psiとストークスベクトルsのパラメータs_{0}，s_{1}，s_{2}の関係は次式で表せられます．

I_{max} = \frac{s_{0} + \sqrt{s_{1}^{2} + s_{2}^{2}}}{2}，I_{min} = \frac{s_{0} - \sqrt{s_{1}^{2} + s_{2}^{2}}}{2}，\psi = \frac{1}{2} \tan^{-1}\frac{s_{2}}{s_{1}}

偏光の計測

偏光の状態を取得するためには，カメラの前に配置した偏光板を回転させながら複数枚撮影した画像（輝度値）からパラメータを推定します．非線形最小二乗法を用いてコサイン関数I(\upsilon)にフィッティングさせることにより求めることができますが計算コストが高いです．一方でストークスベクトルsは線形最小二乗法を用いてパラメータの推定ができるので高速に演算することができます．

偏光板の回転角度を\upsilon_{1},\upsilon_{2}, \cdots, \upsilon_{n}*1と変化させながら観測した輝度値がI(\upsilon_{1}),I(\upsilon_{2}),\cdots,I(\upsilon_{n})となるとき，ストークスベクトルとの関係を行列の形で表すと次式のようになります．

0.5\begin{pmatrix}1 \hspace{1ex} \cos2\upsilon_{1} \hspace{1ex} \sin2\upsilon_{1} \\ 1 \hspace{1ex} \cos2\upsilon_{2} \hspace{1ex} \sin2\upsilon_{2} \\  \vdots \hspace{4ex} \vdots \hspace{5ex} \vdots \hspace{3ex} \\ 1 \hspace{1ex} \cos2\upsilon_{n} \hspace{1ex} \sin2\upsilon_{n}\end{pmatrix} \begin{pmatrix}s_{0} \\ s_{1} \\ s_{2}\end{pmatrix}=\begin{pmatrix}I(\upsilon_{1}) \\ I(\upsilon_{2}) \\ \vdots \\ I(\upsilon_{n})\end{pmatrix}

この式を最小二乗法で解くことでストークスベクトルsのパラメータが求まります．

偏光カメラ

偏光カメラはイメージセンサーの前に0^\circ, 45^\circ, 90^\circ, 135^\circの偏光子が規則的に配置されており，ワンショットで4方向の偏光情報を撮影できます．Sonyは2018年から偏光イメージセンサーの量産を始めました．そのため，最近では偏光カメラが安価に手に入るようになってきています．

www.sony-semicon.co.jp

偏光の応用

偏光カメラ（FLIRのBFS-U3-51S5P-C，SonyのIMX250MZRセンサーを搭載）を用いて撮影を行いました．

偏光度

黒いカメラを撮影しました．通常の撮影では暗くてわかりにくいですが，偏光度\rhoに注目するとカメラの領域だけがはっきりと見えます．影で暗くなっている部分とも見分けがついていることもわかります．

f:id:raccoon15:20191230024249p:plain
f:id:raccoon15:20191230024255p:plain
（左）元画像，（右）偏光度(DoLP)画像

偏光方向（面法線）

アボカドを撮影しました．位相角(0^\circ~180^\circ)の値をHSVで擬似カラー化して表しています．ランダム偏光の光が物体表面で反射すると，フレネル反射に従って反射光は偏光します．偏光した光の位相角や偏光度は光の入射角（＝反射角）と関係があります．そのため，位相角や偏光度から反射面の法線方向を推定することができます（ただし曖昧性の問題があります）．

f:id:raccoon15:20200102002922p:plain
f:id:raccoon15:20200102002926p:plain
（左）元画像，（右）位相角(AoLP)画像

応力・ひずみ

透明のプラスチックのスプーンを液晶ディスプレイ（面発光する偏光光源）の上に置いて，真上から撮影しています．通常では透明なスプーンですが，わずかな歪みが偏光に影響を与えています．

f:id:raccoon15:20191230232634p:plain
f:id:raccoon15:20191230232640p:plain
（左）元画像，（右）位相角+偏光度(AoLP+DoLP)画像

反射除去

池にいる鯉を撮影しました．水面の反射によって外の風景が映り込んでおり，鯉がはっきり見えない場所があります．I_{min}成分を取り出すと，水面の反射成分を除去できています．

f:id:raccoon15:20200103180050p:plain
f:id:raccoon15:20200103180053p:plain
（左）元画像，（右）反射除去(Imin)

拡散反射と鏡面反射成分の分離

シーンに対して完全直線偏光の光で照らしています．拡散反射成分はランダム偏光になるのに対して，鏡面反射成分は偏光が保たれます．そのため，2I_{min}は拡散反射成分，I_{max}-I_{min}は鏡面反射成分に対応しています．

f:id:raccoon15:20200104182500p:plain
f:id:raccoon15:20200104182506p:plain
f:id:raccoon15:20200104182509p:plain
（左）元画像，（中）拡散反射成分(2Imin)，（右）鏡面反射成分(Imax-Imin)

プログラム

偏光カメラ(FLIR, BFS-U3-51S5P-C)で撮影した画像のデモザイキングやストークスベクトルの推定などを行うプログラムをGitHubに置いておきます．この記事に掲載している画像も一緒に置いています．

github.com

参考資料

*1:パラメータの数は3つなので最低3回の撮影が必要です(n>=3)．また，偏光板の回転角度が同じだと不良設定問題になる場合があります．

ダイポールモデルによる半透明物体の散乱特性の計測

質感 ライトトランスポート コンピュータグラフィックス

概要

半透明物体では，光が物体内部に入り込んで反射する表面下散乱という現象が発生します．この現象をダイポールモデルで近似することによって，半透明物体の散乱特性を計測することができます．

今回はろうの散乱特性を計測しました．さらに，計測した結果を用いてろうの質感をCGで再現しました．

半透明物体

私たちの身の周りには半透明物体で溢れています．大理石・皮膚・ろうそく・果物・肉・牛乳など，程度の差はありますが，これら全ては半透明物体としての性質を持っています．

f:id:raccoon15:20191019002914j:plain

ろうにレーザ光を照射した様子．照射された光が物体内部で散乱していることが確認できる．

半透明物体では光が物体内部に入って，別の点から出てくる現象が発生します．このような現象を表面下散乱(Subsurface scattering)と呼びます．表面下散乱では，ある物体表面上の点x_{i}に入った光が物体内で1回ないし複数回反射して別の点x_{o}から出てきます．

f:id:raccoon15:20191019193322p:plain
f:id:raccoon15:20191019193325p:plain
（左）物体表面での反射，（右）物体内部での表面下散乱

表面下散乱が発生する物体の表面上での光の関係を表現した関数をBSSRDF (Bidirectional Scattering Surface Reflectance Distribution Function, 双方向散乱面反射率分布関数)と呼びます．BSSRDFは入射光と出射光の位置x_{i}，x_{o}と方向\vec{\omega_{i}}，\vec{\omega_{o}}の8次元の関数S \left( x_{i},\vec{\omega_{i}}, x_{o}, \vec{\omega_{o}} \right)です．

BSSRDFが既知であれば，その物体をCGで表現できるようになります．しかし，実際に計測しようとすると非常に時間がかかります．また，計測できたとしても莫大なデータ量になります*1．

ダイポールモデル

BSSRDFを表すパラメトリックモデルとして，ダイポールモデル(Dipole model)*2があります．このモデルは，多重散乱*3 *4で，無限平面で，十分に厚みのある*5，均一散乱物体のBSSRDFを良く近似します．

ダイポールモデルではBSSRDFを次式のように構成しています． 

S \left( x_{i},\vec{\omega_{i}}, x_{o}, \vec{\omega_{o}} \right) = \frac{1}{\pi} F_{t}(\eta,\vec{\omega_{i}}) R_{d}(x_{i}, x_{o}) F_{t}(\eta,\vec{\omega_{o}})

\etaは相対屈折率，F_{t}はフレネル透過率です．r=||x_{i}-x_{o}||とすると，散乱項R_{d}は次式で表せられます． 

R_{d}(r)=\frac{\alpha}{4 \pi}\left\{ z_{r}\left(\sigma_{tr} + \frac{1}{d_{r}}\right)\frac{e^{-\sigma_{tr} d_{r}}}{d_{r}^{2}} + z_{v}\left(\sigma_{tr}+\frac{1}{d_{v}}\right)\frac{e^{-\sigma_{tr} d_{v}}}{d_{v}^{2}} \right\}

このとき各変数は次式のようになります． 

d_{r}=\sqrt{r^{2}+z_{r}^2}，d_{v}=\sqrt{r^{2}+z_{v}^2}，z_{r}=\frac{1}{\sigma_{t}^{\prime}}，z_{v}=z_{r}\left(1+\frac{4}{3}A \right)，A=\frac{1+F_{dr}}{1-F_{dr}}，F_{dr}=-\frac{1.440}{\eta^{2}}+\frac{0.710}{\eta}+0.668+0.0636\eta，\sigma_{tr}=\sqrt{3 \sigma_{a} \sigma_{t}^{\prime}}，\sigma_{t}^{\prime}=\sigma_{s}^{\prime}+\sigma_{a}，\sigma_{s}^{\prime}=\sigma_{s}(1-g)，\alpha=\frac{\sigma_{s}^{\prime}}{\sigma_{t}^{\prime}}

この中で，散乱係数(Scattering cofficient)\sigma_{s}，吸収係数(Absorption coefficient)\sigma_{a}，相対屈折率\eta，散乱方向の平均コサイン(Mean cosine of the scattering angle)g*6の4つがパラメータになります． 

例として，皮膚（Skin1[Jensen et al, 2001]）のパラメータを元にR_{d}とrの関係をプロットしました．光が入射された位置から離れるほどR_{d}は指数関数的に減衰します．

f:id:raccoon15:20191019204121p:plain

散乱項R_{d}と距離rの関係（Skin1）

実験

実際に計測を行いダイポールモデルのパラメータを推定します．対象の半透明物体は，ろう（一辺の長さは6.9cm）です．光源として赤色レーザを使いました．カメラはCanon EOS 6DでレンズはSMC TAKUMAR 90mm F2.8です．レーザ光が対象物体に対して垂直に入射するように設置しました．表面下散乱は距離に応じて指数関数的に減衰するので輝度差が大きくなります．そこでHDR撮影を行いました．

f:id:raccoon15:20191021122024j:plain

実験のセットアップ

まず，半透明物体を撮影した画像から，レーザ光の入射位置を原点としたときの距離r_{i}に対応する画素値p_{i}を抜き出しました．レーザ光の放射束\Phiと未知の定数Kを用いるとp_{i}は次式で表せられます．

p_{i} = K \Phi R_{d}(r_{i})

次に，Kと\Phiを求めるために拡散反射板を撮影します．理想的には標準拡散反射板（反射率>0.99）を使いたいのですが，持ってなかった*7のでコピー用紙で代替しました．拡散反射板を撮影した画像の画素値の合計を求めるとK\Phi/Aになります．Aは1画素あたりで撮影された拡散反射板の領域の面積です．コピー用紙に既知の長さの線を印刷しておいたので，それを元にAを求めました．Aが求まったのでK\Phiを求めることができます．

最後にp_{i}/(K\Phi)を計算するとR_{d}が求まります．r_{i}はろうの大きさを定規で測定して，画像中の位置との関係から求めました．

f:id:raccoon15:20191020013827j:plain
f:id:raccoon15:20191020013831j:plain
撮影の対象．（左）拡散反射板（コピー用紙で代替），（右）半透明物体（ろう）．

結果

レーザ光の入射位置を原点としたときの距離rを横軸に，計測した散乱項R_{d}を縦軸にしてプロットしました．原点から離れるとR_{d}は小さくなっていることが確認できます．原点に近い点r<1では，R_{d}が極端に大きくなっています．この点では物体表面での直接反射成分が含まれているためだと考えられます．

f:id:raccoon15:20191020151457p:plain

計測したろうの散乱項R_{d}と距離r．推定したパラメータを用いて曲線を描画した．

ダイポールモデルに基づいてパラメータの推定を行います．ここでは等方散乱を仮定するのでg=0とします．相対屈折率は多くの物体で\eta=1.3となることが知られているのでこの値を用います．非線形最適化でパラメータを推定すると，\sigma_{s}^{\prime}=0.112/mm，\sigma_{a}=0.00711/mmとなりました．

推定したパラメータを用いてレンダリングを行いました．レンダラはMitsubaを使いました．

f:id:raccoon15:20191019015226p:plain
f:id:raccoon15:20191019015230p:plain
レンダリング結果．（左）拡散反射，（右）ダイポールモデル（計測したろう）．

表面下散乱の生じない拡散反射と比べると，ダイポールモデルでは陰影が少なくなり，半透明物体の質感が表現できていると思います． 

立方体以外の複雑なシーンにも適用してみました．

f:id:raccoon15:20191019114738p:plain
f:id:raccoon15:20191020163229p:plain
ろうの質感を表現できる．

参考資料

*1:100x100の画素数で1度単位の計測を行うと，入射光位置は100x100通り，出射光位置は100x100通り，入射光方向は360x90通り，出射光方向は360x90通り，となる．全ての組み合わせは10^{16}通り．8bitで記録しても100PB（ペタバイト）必要になる．

*2:H. W. Jensen, S. R. Marschner, M. Levoy, and P. Hanrahan, “A Practical Model for Subsurface Light Transport”, Proc. SIGGRAPH2001, 2001.

*3:物質内で何度（2回以上）も反射する散乱を多重散乱．1回だけの反射は単一散乱と呼ぶ．

*4:単一散乱も近似できるモデルとして，Directional dipole modelが提案されている．J. R. Frisvad, T. Hachisuka, and T. K. Kjeldsen, “Directional dipole model for subsurface scattering,” ACM Trans. Graph., 2014.

*5:特定の厚みを持つ物体を近似できるモデルとして，Multipole modelが提案されている．C. Donner H. W. Jensen, “Light Diffusion in Multi-Layered Translucent Materials”, Proc. SIGGRAPH2005, 2005.

*6:散乱の方向を表すフェーズ関数のパラメータ．-1.0~1.0の範囲の値をとる．g=0のとき等方散乱，g>0のとき前方散乱，g<0のとき後方散乱になる．

*7:10万円くらいするみたいですね．Labsphere Spectralon, Diffuse Reflectance．

レンジファインダによる3次元計測

3次元計測 ライトトランスポート

概要

カメラとプロジェクタのペアを使って，物体の3次元計測を行いました．プロジェクタから投影するパターンはグレイコードを使いました．

レンジファインダ

3次元計測を行う手法としては2台のカメラを用いたステレオ法が有名です．2台のカメラの視差を利用することで，三角測量の原理に基づき距離を算出することができます．この手法では，テクスチャのない部分では対応点付けが困難なため，距離計測が不安定になります．

レンジファインダは，カメラとプロジェクタのペアで距離を計測できます．プロジェクタから光のパターンを投影して，カメラで撮影します．すると，テクスチャのない物体でも，対応点付け問題が簡単になります． 

f:id:raccoon15:20190503210716p:plain

パターン投影することで距離を計測できる（グレイコードの例）

3次元計測の手順

カメラとプロジェクタの対応関係の獲得

3次元計測を行うためには，プロジェクタから出力された光によって光っている物体上の点が，カメラの画像上のどの座標で撮影されるかという対応関係を求める必要があります．

もっとも簡単な方法は，プロジェクタからある1点だけを光らせ，カメラで撮影することです．これを全ての点に対して行うことで，プロジェクタとカメラの対応関係が求まります（スポット光投影法）．しかし，この方法では非常に時間がかかります．

計測時間を短縮するための手法として，グレイコードパターンを投影する方法があります．この手法では白黒のパターンを複数枚投影して撮影し，デコードを行うことで，カメラとプロジェクタの対応関係を求めることができます．例えば，横幅が1024ピクセルのプロジェクタであれば，最低10枚（反転パターンを含めると計20枚）のパターンを投影して撮影することで対応関係（カメラの2次元座標とプロジェクタの1次元座標との対応）を求めることができます．

blog.k-tanaka.me

カメラとプロジェクタのキャリブレーション

カメラのキャリブレーションをすると，次式で表せられるカメラパラメータCを獲得できます．(X, Y, Z)は世界座標系でのある点で，(X_{c}, Y_{c})はその点に対応するカメラ座標系での画像の二次元座標です．

h \begin{pmatrix} X_{c} \\ Y_{c} \\ 1 \end{pmatrix} = \begin{pmatrix} C_{11}  C_{12}  C_{13}  C_{14} \\ C_{21}  C_{22}  C_{23}  C_{24} \\ C_{31} C_{32} C_{33} C_{34} \end{pmatrix}\begin{pmatrix} X \\ Y \\ Z \\ 1 \end{pmatrix}

プロジェクタのキャリブレーションをすると，次式で表されるプロジェクタパラメータPを獲得できます．(X, Y, Z)は世界座標系でのある点で，X_{p}はその点に対応するプロジェクタ座標系でのx軸方向の座標です．

h \begin{pmatrix} X_{p} \\ 1 \end{pmatrix} = \begin{pmatrix} P_{11}  P_{12}  P_{13}  P_{14} \\ P_{21}  P_{22}  P_{23}  P_{24}\end{pmatrix}\begin{pmatrix} X \\ Y \\ Z \\ 1 \end{pmatrix}

raccoon15.hatenablog.com

距離の計算 

カメラパラメータCとプロジェクタパラメータPを用いて，行列F，Q，Vを定義します．

F = \begin{pmatrix}C_{34}X_{c}-C_{14} \\ C_{34}Y_{c}-C_{24} \\ P_{24}X_{p}-P_{14}\end{pmatrix}

Q = \begin{pmatrix}C_{11}-C_{31}X_{c} \ C_{12}-C_{32}X_{c} \ C_{13}-C_{33}X_{c} \\ C_{21}-C_{31}Y_{c} \ C_{22}-C_{32}Y_{c} \ C_{23}-C_{33}Y_{c} \\ P_{11}-P_{21}X_{p} \ P_{12}-P_{22}X_{p} \ P_{13}-P_{23}X_{p}\end{pmatrix}

V = \begin{pmatrix}X \\ Y \\ Z\end{pmatrix} 

カメラ座標(X_{c}, Y_{c})とプロジェクタ座標X_{p}に対応する物体座標(X, Y, Z)は次式で計算できます．

V=Q^{-1}F

この計算を画像上の全ての座標で行うことで，物体の3次元座標を求めることができます．

実験結果

グレイコードパターンを投影することで，プロジェクタとカメラの対応点を求め，物体の3次元計測を行いました．グレイコードは，ポジのパターンを11枚とネガのパターンを11枚投影しました．

f:id:raccoon15:20190509185411j:plain

実験の様子

対象物体は，紙コップとキャリブレーションボックスです．

f:id:raccoon15:20190509185446j:plain
f:id:raccoon15:20190509185451j:plain
対象物体

3次元計測を行なった結果です．

f:id:raccoon15:20190509190303p:plain
f:id:raccoon15:20190509190445p:plain
計測結果


Active Stereo Result

基準立方体を用いたカメラキャリブレーション

キャリブレーション

概要

立方体を用いてカメラのキャリブレーションを行い，カメラパラメータを求めました．立方体を用いることでワンショットでキャリブレーションができます．

カメラキャリブレーション

幾何学的カメラキャリブレーションの手法としてはZhangのチェスボードを用いる方法が有名です．この手法ではチェスボードを動かしながら複数枚撮影することで，キャリブレーションをすることができます．チェスボードは手軽に用意することができます．また，OpenCVやMATLABでは専用のプログラムが実装されています．

チェスボードを用いない手法としては，形が既知の物体を用いる方法があります．

f:id:raccoon15:20190428230200j:plain

基準立方体

今回は上の図のような基準立方体を使います．立方体の大きさは8×8×8cmです．チェスボードと比較して，基準立方体を使うメリットは2つあります．

1. ワンショットでキャリブレーションを行うことができる

基準立方体をカメラで撮影すると，3面を同時に撮影することができます．そのため，チェスボードのように動かしながら複数枚撮影する必要がありません．

2. プロジェクタのキャリブレーションにも適用しやすい

プロジェクタのキャリブレーションを行う場合は，パターンを投影してカメラで撮影します．このとき黒い領域では，投影したパターンが暗くなるので観測しにくくなります．この基準立方体では，各面には辺に沿って黒いマーカをつけています．このマーカを対辺どうしで直線で結んだ時の縦線，横線が交わる点を基準点として獲得できます．このとき基準点は全て白い領域なので，投影したパターンを観測しやすいという利点があります．

立方体の世界座標

まず，基準立方体の座標系について設定します．立方体のある頂点Oを原点として定め，世界座標系で(0, 0, 0)としました．辺に沿って，x軸，y軸，z軸を設定しました．立方体の大きさは既知なので，各面上の点の座標がわかります．

f:id:raccoon15:20190503014319j:plain

基準立方体の世界座標

立方体のカメラ座標の獲得

基準立方体の世界座標と対応するカメラ座標を求めます．

基準立方体を撮影した画像を用意します．背景は模様が無いようにすると，後の処理が楽になります．

f:id:raccoon15:20190502234545p:plain

撮影した基準立方体

撮影した画像を2値化して，黒いマーカがブロブとして検出できるようにします．

f:id:raccoon15:20190502234558p:plain

2値化の結果

ラベリングを行い，黒いマーカの位置を求めます．

f:id:raccoon15:20190502234603p:plain

黒いマーカの検出

対辺どうしのマーカを直線で結んだとした時の縦線，横線が交わる格子点を基準点とします．この立方体では，一辺に7個のマーカがあるので，基準点は一面に49点．それが三面あるので合計で147点数の基準点が得られます．

f:id:raccoon15:20190502234610p:plain

格子点を求める

カメラパラメータの計算

世界座標系でのある点(X, Y, Z)と，カメラ座標系での画像の二次元座標(X_{c}, Y_{c})との関係は次式で表すことができます．

h \begin{pmatrix} X_{c} \\ Y_{c} \\ 1 \end{pmatrix} = \begin{pmatrix} C_{11}  C_{12}  C_{13}  C_{14} \\ C_{21}  C_{22}  C_{23}  C_{24} \\ C_{31} C_{32} C_{33} C_{34} \end{pmatrix}\begin{pmatrix} X \\ Y \\ Z \\ 1 \end{pmatrix}

この式を展開して，媒介変数hを消去すると，

(C_{11}-C_{31}X_{c})X+(C_{12}-C_{32}X_{c})Y+(C_{13}-C_{33}X_{c})Z = C_{34}X_{c}-C_{14}

(C_{21}-C_{31}Y_{c})X+(C_{22}-C_{32}Y_{c})Y+(C_{23}-C_{33}Y_{c})Z = C_{34}Y_{c}-C_{24}

式より，(X, Y, Z)と(X_{c}, Y_{c})の対応が既知の点が1点あれば2つの方程式を作ることができます．求めたいカメラパラメータCは3×4の行列なので，未知数は12個です．そのため最低6組の点があれば方程式を12個作ることができ，カメラパラメータCを求めることができます．

実際には6組以上の点を使い，最小二乗法を適用することでパラメータを求めるます．

プロジェクタキャリブレーション

カメラパラメータを求める式を用いることで，プロジェクタのパラメータを求めることができます．このとき必要になるのは，プロジェクタ座標系のある点の座標と，それに対応する世界座標系の点の座標です．プロジェクタは光を出力する装置なので，カメラのように直接撮影して座標を求めることができません．そのため，プロジェクタから投影したパターンをカメラで撮影することで座標の対応関係を獲得する方法がよく使われます．もっとも簡単な方法は，プロジェクタからある1点だけを光らせ，それをカメラで撮影することです．これを全ての点に対して行うことで，プロジェクタとカメラの対応関係を求まります．しかし，この方法では非常に時間がかかるので，グレイコードパターンを投影する方法などが用いられます．

The Imaging Source社のカメラをPythonから動かす

機材

概要

The Imaging Source社のマシンビジョン用のカメラをPythonから動かせたのでメモしておきます． 

f:id:raccoon15:20190420234250j:plain

The Imaging Source社のカメラからキャプチャを行なった

実行環境

やり方

1.ドライバのインストール

カメラのドライバをインストールします．

https://www.argocorp.com/software/DL/tis/driver.html

2.コントローラのダウンロード

The Imaging Source社のGitHubからコントローラをダウンロードします．

github.com

git clone https://github.com/TheImagingSource/IC-Imaging-Control-Samples.git

Pythonのサンプルがあるディレクトリへ移動します．

cd '.¥IC-Imaging-Control-Samples¥Python¥Open Camera, Grab Image to OpenCV¥'

サンプルプログラムを実行します．

python .¥tis-OpenCV.py

実行すると，カメラからキャプチャされた画像が表示されます．

 

プログラム

自分なりに改造した撮影用プログラムをGitHubに置いておきます．

このプログラムでは，露光時間・ゲインの設定，複数枚撮影によるノイズ低減，簡易のHDR画像撮影機能を備えた自作関数を実装しています．

github.com