チンダル現象を用いた実物体の三次元形状及 び双方向

修士論文
チンダル現象を用いた実物体の三次元形状及
び双方向反射率分布関数推定
指導教員
美濃 導彦 教授
京都大学大学院情報学研究科
修士課程知能情報学専攻
小山 亮
平成 23 年 2 月 8 日
i
チンダル現象を用いた実物体の三次元形状及び双方向反射率分布
関数推定
小山 亮
内容梗概
近年，映画，アニメーションなどで実物体の見えを忠実に再現した写実的な
CG が頻繁に用いられるようになった．しかし，写実的な CG の生成には多大
な人的コストがかかるため，これを自動で生成するという要求が高まってきて
おり、実物体の写実的な CG を自動で生成するのに必要な三次元形状及び表面
反射特性を実物体より計測することが求められている．そのうち表面反射特性
については一般に双方向反射率分布関数 (BRDF) を求めることで推定される．
従来より物体表面に光を照射しその反射光をカメラで観測した画像から物体
の三次元形状及び表面反射特性を推定する手法が多く提案されている．しかし，
このような手法全般に共通する問題点として，反射光がカメラで観測できない
黒色物体や金属物体の場合は三次元形状を推定することができないという問題
がある．
このような従来研究に対して，我々はこれまで反射光がカメラで観測できな
いために計測が困難であった物体についても計測が可能となる手法を提案する．
提案手法では，この問題に対処するために，チンダル現象という光学現象に
着目する．チンダル現象が生じない場合では光の道筋は観測できないが，多数
の微粒子で構成される媒体中ではチンダル現象による散乱が発生し，入射光及
び反射光が媒体中を通過する様子が観測できる．
これにより，通常の環境では観測ができなかった物体表面での反射光につい
てもカメラで観測できるようになる．入射光が媒体中を通過する様子を観測す
ることによって光の反射位置，つまり物体表面の位置を知ることができるため，
黒色物体のように光をほとんど反射しないような物体に対しても形状計測が可
能となる．さらに，物体表面で反射したあらゆる方向での反射光を観測するこ
とで数式モデルを用いない BRDF を推定する．
本研究では BRDF を Lambert モデルと Torrance-Sparrow モデルによって表
現し，そのパラメータを獲得する．金属物体と黒色物体を計測対象とした実験
により提案手法の有効性を確認した．
ii
The Estimation of the Shape and Bidirectional
Reflectance Distribution Function of Objects under
Tyndall Effect
Akira KOYAMA
Abstract
Recently, realistic CG is used for movies, news shows, and animations. Creating realistic CG by hand costs enormously, so there are increasing demands
for creating realistic CG automatically. 3D shape and reflection characteristics
measured from real objects are necessary for creating realistic CG. The reflection
characteristics are usually described as Bidirectional Reflectance Distribution
Function (BRDF).
Some methods for estimating the 3D shape and BRDF of objects are proposed. They estimated the shape and BRDF by irradiating the object and
capturing the reflected light with cameras. However, they cannot estimate the
3D shape of the objects such as black objects or metal objects if cameras cannot
capture the reflected light. We propose the method to estimate the 3D shape
of such objects.
To overcome the problem of no capturing the reflected lights by the cameras,
we take advantage of Tyndall effect. Under Tyndall effect, we can observe
directly the incident and reflected light from any view points. This is because
the particles let the light scatter.
Observing the line of incident light under Tyndall effect, we can estimate
the 3D shape of objects as the reflection point where the light broke off. In
addition, we can capture scattered light of the reflection at the surface of objects and estimate the BRDF of objects as non-parametric BRDF model. In
this study, we describe BRDF by Lambert model and Torrance-Sparrow model
and estimate the parameters. We evaluated the effectiveness of the proposed
method, estimating the 3D shape and reflection characteristics of a metal object
and a black object.
チンダル現象を用いた実物体の三次元形状及び双方向反射率分布
関数推定
目次
第1章
はじめに
1
第2章
関連研究
5
2.1
三次元形状推定に関する関連研究 . . . . . . . . . . . . . . . . . . . . .
5
2.2
BRDF 推定に関する従来研究 . . . . . . . . . . . . . . . . . . . . . . . .
6
2.3
三次元形状と表面反射特性の同時推定に関する関連研究 . . . . .
9
2.4
関連研究における問題点に対する解決策 . . . . . . . . . . . . . . . . . 10
第3章
チンダル現象下での三次元形状推定
11
3.1
観測方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2
カメラキャリブレーション . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.3
第4章
3.2.1
カメラキャリブレーションの概要 . . . . . . . . . . . . . . . . 11
3.2.2
内部キャリブレーション . . . . . . . . . . . . . . . . . . . . . . 13
3.2.3
外部キャリブレーション . . . . . . . . . . . . . . . . . . . . . . 14
三次元形状推定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
チンダル現象下での BRDF 推定
20
4.1
チンダル現象下での BRDF 推定手法の概要 . . . . . . . . . . . . . . 20
4.2
散乱現象 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.3
シミュレート画像生成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.3.1
光線追跡 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.3.2
単位空間の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.3.3
シミュレート画像生成アルゴリズム . . . . . . . . . . . . . . 25
4.4
評価関数の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.5
パラメータ探索による BRDF 推定 . . . . . . . . . . . . . . . . . . . . . 29
第5章
実験
31
5.1
計測装置 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.2
誤差評価関数の正当性の評価 . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.2.1
パラメータ 1 種類を変化させた場合 . . . . . . . . . . . . . . 32
5.2.2
5.3
パラメータ 2 種類を変化させた場合 . . . . . . . . . . . . . . 33
観測画像を用いた三次元形状，BRDF パラメータ，散乱パラメー
タ推定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
第6章
5.3.1
黒色物体 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.3.2
金属物体 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
おわりに
謝辞
53
54
参考文献
55
第1章
はじめに
近年映画，ニュース番組，アニメーションなどで実物体の見えを忠実に再現
した写実的な CG が頻繁に用いられるようになった．しかし，写実的な CG の
生成には多大な人的コストがかかるため，これを自動で生成するという要求が
高まってきている．実物体の写実的な CG を自動で生成するにはその物体の三
次元形状及び表面反射特性が必要であり，これを実物体より計測することが求
められている．
実物体の三次元形状を推定する手法については，従来より物体表面に光を照
射しその反射光をカメラで観測した画像から光の反射位置を推定する手法が多
く提案されている [1][2][3]．しかし，このような手法全般に共通する問題点とし
て，反射光がカメラで観測できない物体の場合は三次元形状を推定することが
できないという問題がある．例えば，反射光が一部の方向ではほとんど観測で
きない金属物体や反射光がほとんど発生しない黒色物体には適用できない．
一方，表面反射特性については一般に双方向反射率分布関数 (Bidirectional
Reflectance Distribution Function, BRDF) で表現されることが多い．BRDF と
は，入射方向及び反射方向に対する反射率 (=反射光の光量/入射光の光量) を記
述する関数である．
BRDF を推定している従来手法としては，物体表面に光を照射し，様々な入
射方向及び反射方向での反射率を計測して様々な物体の BRDF を推定する研究
が提案されている [4][5][6]．しかしながら町田らの手法 [4] では物体の三次元形
状が既知である必要があり，Li らの手法 [5] 及び向川らの手法 [6] では物体表面
上のある一点の BRDF を推定することを目的としているため，計測対象が平面
である必要があるなど計測装置の物理的制約から計測できる物体形状に大きな
制限がある．
また，三次元形状と BRDF を同時に計測する手法がいくつか提案されている
[7] ものの，三次元形状推定における問題点である反射光がカメラで観測できな
い場合についての対応は，この従来研究でも解決できていない．
これらの従来研究に対して我々は，これまで反射光がカメラで観測できないた
めに計測が困難であった金属物体や黒色物体についても計測が可能となる手法
を提案する．また，物体の三次元形状だけでなく数式モデルを用いない BRDF
も同時に獲得する．
1
提案手法では，反射光がカメラで観測できないという問題に対処するために，
チンダル現象という光学現象に着目する．チンダル現象とは光が多数の微粒子
で構成される媒体中を通過する際，その微粒子で光が散乱し光の道筋が見える
ようになる現象である．
チンダル現象の例として，図 1，図 2 にそれぞれ水に洗剤を混ぜた混合物，暗
室で密閉し煙を充満させた空間にレーザを照射した様子を示す．チンダル現象
が生じない場合では光の道筋は観測できないが，洗剤や煙といった多数の微粒
子で構成される媒体中ではチンダル現象による散乱が発生し，入射光及び反射
光が媒体中を通過する様子が観測できる．
これにより，通常の環境ではカメラに入射しないため観測ができなかった物
体表面での反射光についてもカメラで観測できるようになる．また，入射光が
媒体中を通過する様子を観測することによって光の反射位置，つまり物体表面
の位置を知ることができるため，黒色物体のように光をほとんど反射しないよ
うな物体に対しても形状計測が可能となる．さらに，物体表面で反射したあら
ゆる方向での反射光がカメラで一度に観測されるため，本研究ではこれを利用
して数式モデルを用いない BRDF を推定する．
図 1: 水と洗剤の混合物に対するレー 図 2: 煙が充満した空間でのレーザ
ザ照射
照射
本研究では，計測対象となる物体を散乱媒体中に設置し，物体表面にレーザ
光を照射した際に入射光及び反射光がチンダル現象によって散乱する様子を物
体を取り囲むように設置した複数台のカメラで観測する．そのようにして観測
された画像の例を図 3 に示す．このような画像からレーザ光が物体表面に到達
し，様々な方向に反射する様子を観測する．
2
三次元形状推定に関して，まず観測画像から画像処理を用いて入射方向，反
射位置を推定し，各カメラの画像での推定結果からステレオ視により三次元空
間での入射方向，反射位置を求める．この反射位置が物体表面位置の一点とな
るので，チンダル現象下で物体の様々な表面位置にレーザを当てることで物体
表面位置を表す点群，すなわち三次元形状が推定できる．
(a)
(b)
(c)
図 3: 観測画像
三次元形状を推定する際は画像中でどこが明るくてどこが暗いかといった輝
度値の大小の情報があれば十分であったが，BRDF を推定する際は輝度値の大
小だけでなく，入射光に対する反射光の光量の比を画像の輝度値から推定する
必要がある．画像で観測される輝度値は，カメラに直接入る反射光の他に，チ
ンダル現象による散乱光の影響を受けており，輝度値から直接 BRDF を推定す
るためには，BRDF と同時に散乱現象の影響も推定しなければならない．
散乱現象に関して，その振る舞いは散乱パラメータによって表される．本研
究ではその散乱パラメータ及び BRDF を同時に推定する．これらを本研究では
パラメータ探索によって推定する．具体的には，推定した光源と，与えられた
散乱パラメータ及び各反射方向への反射率の下で各カメラで観測される画像を
光線追跡法を用いて生成する (以降，これをシミュレート画像と呼ぶ)．誤差評
価関数によりシミュレート画像と元の観測画像の誤差を計算し，最も評価関数
の値が小さくなるような BRDF 及び散乱パラメータを推定する．提案手法では
BRDF を Lambert モデルと Torrance-Sparrow モデルによって表現し，これら
のパラメータを推定することで計測対象の BRDF を得る．
本研究では上記の手法を用いてレーザ光を物体に照射したときのレーザ光反射
位置及び BRDF，散乱パラメータを推定する．以下，2 章で三次元形状や BRDF
の推定に関する関連研究を挙げる．3 章では三次元形状推定方法について述べ
3
る．4 章ではシミュレート画像の生成方法及び BRDF パラメータ，散乱パラメー
タを推定する方法について述べる．5 章では提案手法の有効性を検証するため
の実験及びその結果・考察を述べる．6 章で本研究のまとめと今後の課題につ
いて述べる．
4
第2章
関連研究
本研究に関連する研究としては，実物体の三次元形状のみを推定する研究，
BRDF のみを推定する研究，三次元形状と BRDF を同時に推定する研究がある．
本章では，これらの関連研究やその問題点について 2.1-2.3 節で述べ，その問題
点の解決策について 2.4 節で述べる．
2.1
三次元形状推定に関する関連研究
本節では，物体表面に光を照射しその反射光をカメラで観測した画像から推
定した光の反射位置をもとに物体の三次元形状を推定している関連研究につい
て述べる．
その一つとして，レーザレンジファインダの Time-of-Flight 方式距離センサ
を用いた三次元形状計測手法 [2] がある．その計測原理は図 4 に示すように，セ
ンサから照射されたレーザ光が物体表面で反射して再びセンサに戻るまでの時
間を計測することで物体との距離を得るというものである．これを用いること
で物体の三次元形状を推定することができる．
図 4: レーザレンジファインダの計
測原理
一本のレーザ光ではなく，平面状のスリット光を物体に投射しその三次元形
状を推定する光切断法 [3] がある．図 5 のようにスリット光を物体に投射すると，
カメラではそれを一本の明るい縞として観測することができ，その縞の曲がり
具合からその縞の部分に相当する物体表面の形状が得られる．そのスリット光
をずらしながら観測を行うことで，物体の三次元形状を推定することができる．
また，プロジェクタを用いたワンショットスキャン法 [1] と呼ばれる手法があ
る．この手法では，図 6 にあるような 2 色のラインパターン (縦，横) をプロジェ
クタを用いて物体に投射する．すると，図 7 のように対象物体の物体表面上に
5
図 5: 光切断法の計測原理
縦横のラインが投影されるので，その様子をカメラで観測し縦横のラインの交
点を推定することで物体の三次元形状を推定している．この手法はレーザレン
ジファインダを用いた手法と比べて三次元形状を表す物体表面位置の点群の密
度は粗いが，一度の撮影で三次元形状を推定することができる．
図 7: ラインパターンの投射に
図 6: 縦横のラインパターン
よる形状計測システム
これらの手法では照射した光の物体表面での反射光をカメラで観測すること
で三次元形状を推定しているが，このような手法全般的に共通する問題点とし
て，反射光がカメラで観測できない物体の場合は三次元形状を推定することが
できないという問題がある．例えば，反射光がほとんど発生しない黒色物体や
反射光が一部の方向ではほとんど観測できない金属物体には適用できない.
2.2
BRDF 推定に関する従来研究
一方，BRDF を推定している従来手法としては，物体表面に光を照射し，様々
な入射方向及び反射方向での反射率を計測して様々な物体の BRDF を推定する研
6
究が提案されている [4][5][6]．BRDF とは，入射方向 (θi , ϕi ) と反射方向 (θr , ϕr )
に対する反射率 (=反射光の光量/入射光の光量) を記述する関数である (0 ≤
θi <
π
,0
2
≤ θr <
π
,0
2
≤ ϕi < 2π, 0 ≤ ϕr < 2π ，図 8)．このときの反射率
fr (0 ≤ fr ≤ 1) は式 (1) のような入射方向 (θi , ϕi ) と反射方向 (θr , ϕr ) の関数とし
て表すことができる．
fr = f (θi , ϕi , θr , ϕr )
(1)
z
౉኿శ
෻኿శ
θr
θi
φr
y
φi
x
図 8: BRDF を表現する 4 つの角度
パラメータ
町田らの研究 [4] では BRDF を拡散反射成分と鏡面反射成分に分離し，それ
ぞれを Lambert モデル，Torrance-Sparrow モデル [8] で近似した以下の式 (2) で
表し，BRDF を近似した反射モデルのパラメータを推定することで BRDF を推
定している．
ks
fr = kd (N · L) +
exp
N·V
(
α2
− 2
2σ
)
·G·F
(2)
ここで，G, F はそれぞれ式 (3)，(4) で表され，
2(N · H)(N · V) 2(N · H)(N · L)
,
}
V·H
V·H
(3)
[L · N(g + L · N) − 1]2
1 (g − L · N)2
{1
+
}
2 (g + L · N)2
[L · N(g − L · N) + 1]2
(4)
G = min{1,
F =
g は式 (5) で表される．
g=
√
n2 + (L · N)2 − 1
(5)
この式において，拡散反射係数 kd ，鏡面反射係数 ks ，表面粗さ係数 σ ，物体
の屈折率 n の 4 つの BRDF パラメータの値が分かれば，任意の入射方向，反射
7
方向における反射率を求めることができる．また，図 9 の通り，L は入射方向
ベクトル，V は反射方向ベクトル，N は法線方向ベクトル，H は L と V を二等
分する方向のベクトルであり，以上のベクトルは全て単位ベクトルとする．ま
た α は N と H の間の角度である．
なお，町田らの研究 [4] では，この Torrance-Sparrow 反射モデルを簡略化し
た式 (6) のパラメータ (拡散反射係数 kd ，鏡面反射係数 ks ，表面粗さ係数 σ) を
求めていた．
ks
fr = kd (N · L) +
exp
N·V
(
α2
− 2
2σ
)
(6)
図 9: Lambert モデル，Torrance-Sparrow
モデルにおける幾何学的な関係
Li らの研究 [5] では，光源とカメラをそれぞれ物体の全周に動かしながら様々
な入射方向及び反射方向での反射率を計測し BRDF を推定している．向川らの
手法 [6] では，楕円鏡とプロジェクタを用いた装置を用いて BRDF を推定して
いる．図 10 に示すような楕円鏡では，二つある焦点のうち一方の焦点から出た
光が側面で反射した後他方の焦点に向かうという性質があり，その性質を利用
する．図 11 にあるようにまず光学的に焦点位置にあるプロジェクタの一画素の
みを白くして一本の入射光を照射すると，楕円鏡の側面で反射して他方の焦点
にある物体表面に光が到達する．そして物体表面で様々な方向に反射した光は
楕円鏡の側面で反射して全て再び焦点位置 (光学的にはプロジェクタと同じ位
置) にあるカメラに到達するため，一度の観測で様々な反射方向での反射率の分
布を得ることができ，Li らの手法と比べて大幅に計測時間を短縮している．
8
しかしながら町田らの手法 [4] では物体の三次元形状が既知である必要があ
り，Li らの手法 [5] 及び向川らの手法 [6] では物体表面上のある一点の BRDF を
推定することを目的としているため，計測対象が平面である必要があるなど計
測装置の物理的制約から計測できる物体形状に大きな制限がある．
図 10: 楕円鏡の二つある焦点の一方
の焦点から出た光は側面で反射し他 図 11: 楕円鏡とプロジェクタを用い
方の焦点に向かう
2.3
た BRDF 計測装置
三次元形状と表面反射特性の同時推定に関する関連研究
本研究の他にも，三次元形状と BRDF を同時に計測する手法が提案されてい
る．小俣らの研究 [7] では，図 12 のようなターンテーブルに乗った物体に光を
当て，反射光を真上に設置したカメラで観測を行う．物体に照射する平行光源
を回転させながら取得した観測画像から，図 13 のような物体の回転による明度
変化のグラフが得られる．物体表面が完全拡散面である場合はこの曲線が正弦
波曲線となることを利用し，この中から拡散反射成分しか含まないデータを抽
出した．さらに BRDF を数式モデルで近似した Phong の鏡面反射モデル [9] を
用いて法線方向を含む三次元形状及び表面反射特性を推定している．
しかしながら，この手法で BRDF として推定されるのは，Phong の反射モデ
ル [9] のパラメータであり，計測対象としては Phong の反射モデルに合致する物
体に限られている．なお，数式モデルを用いずに BRDF を推定する手法として
2.2 節で挙げた向川らの研究 [6] があるが，2.2 節で述べた通り計測できる物体
形状に大きな制限がある．また，三次元形状推定の問題点である反射光がカメ
ラで観測できない場合についての対応は，この従来研究でも解決できていない．
9
図 12: ターンテーブルを用いた物体 図 13: 物体の回転によるカメラでの
の三次元形状及び BRDF 計測装置 明度変化
2.4
関連研究における問題点に対する解決策
これらの従来研究に対して我々は，これまで反射光がカメラで観測できないた
めに計測が困難であった黒色物体や金属物体についても計測が可能となる手法
を提案する．また，物体の三次元形状だけでなく数式モデルを用いない BRDF
も同時に獲得する．
提案手法では，反射光がカメラで観測できないという問題に対処するために，
チンダル現象という光学現象に着目する．チンダル現象とは光が多数の微粒子
で構成される媒体中を通過する際，その微粒子で光が散乱し光の道筋が見える
ようになる現象である．チンダル現象が生じない場合では光の道筋は観測でき
ないが，多数の微粒子で構成される媒体中ではチンダル現象による散乱が発生
し，入射光及び反射光が媒体中を通過する様子が観測できる．
これにより，通常の環境ではカメラに入射しないため観測ができなかった物
体表面での反射光についてもカメラで観測できるようになる．また，入射光が
媒体中を通過する様子を観測することによって光の反射位置，つまり物体表面
の位置を知ることができるため，黒色物体のように光を反射しないような物体
に対しても形状計測が可能となる．さらに，物体表面で反射したあらゆる方向
での反射光がカメラで一度に観測されるため，本研究ではこれを利用して数式
モデルを用いない BRDF を推定する．
10
第3章
チンダル現象下での三次元形状推定
本章では本研究における実物体の三次元形状の推定手法について述べる．ま
ず提案手法で想定している計測対象の観測方法に関して 3.1 節で述べる．提案
手法では複数台のカメラを用いて観測を行う．次にカメラ間の位置関係や各カ
メラの焦点距離などをキャリブレーションする方法に関して 3.2 節で述べる．最
後に三次元形状推定のアルゴリズムに関して 3.3 節で述べる．
3.1
観測方法
本節では本研究での観測方法について述べる．計測環境としては，チンダル
現象が発生するような散乱媒体で満たされている空間を考える．この際，散乱
媒体の濃度が不均一な場合，後述する提案手法において観測画像から推定すべ
き三次元形状や BRDF を正しく推定できないということから，散乱媒体の密度
は空間内で均一であるとする．その空間において，静止した計測対象の物体表
面上の一点にレーザ光を照射する．このとき，レーザー光とレーザ光の反射光，
およびレーザ光と反射光がチンダル現象によって散乱した光以外の光は観測さ
れないものとし，対象物体としてはレーザ光が常に物体表面上の一点のみで反
射し物体表面上で多重反射することがないような凸物体を用いる．その様子を
K 台のカメラで複数方向から観測する．カメラの内部パラメータ・外部パラメー
タについては 3.2 節で後述する方法により予め計測しておく．
なお，現実には一度散乱した光が媒体中で再び散乱する多重散乱現象が発生
しているが，一次散乱光に比べて多重散乱光の光量が少ないため，画像中で多
重散乱がほとんど観測されない程度に散乱媒体の濃度やカメラのシャッター開
放時間の調整を行うことによって多重散乱光を無視し一次散乱光のみを扱う．
カメラで観測した画像について，以下ではカメラ k での観測画像を Ikobs (k =
1, · · · , K) とし，その画素 (u, v) での輝度値を ikobs (u, v) とする (u = 1, 2, · · · , Wobs ; v =
1, 2, · · · , Hobs ; 観測画像の解像度は Wobs × Hobs )．
3.2
3.2.1
カメラキャリブレーション
カメラキャリブレーションの概要
カメラを用いて物体を撮影する際，そのカメラの焦点距離やレンズの歪み及
びカメラの位置や向きなどの要因によりカメラ内での物体の位置，向き，大き
11
図 14: 計測装置
さなどが変化する．複数台のカメラを用いて三次元計測を行う際には予め各カ
メラの焦点距離やレンズの歪み，カメラ間の位置関係を計算しておく必要があ
る．それらの情報を数値で表したものがカメラパラメータであり，図 16 のチェ
スパターンなどサイズが既知の物体を撮影しカメラパラメータを推定すること
をカメラキャリブレーションと呼ぶ．
カメラ座標系のモデルとしては図 15 に示すピンホールカメラモデルを用いる．
XY 平面がスクリーンに平行となるように X 軸，Y 軸をとり，X 軸，Y 軸，Z
軸がこの順番で右手系となるように Z 軸をとったとき，カメラの焦点 O を原点
とした O − XY Z 座標系をカメラ座標系と呼ぶ．Z 軸はカメラ光軸とも呼ばれ
る．また，Z 軸とスクリーンとの交点 c を画像中心といい，スクリーンが uv 平
面となり，u 軸，v 軸，Z 軸がこの順番で右手系となるように u 軸，v 軸をとっ
たとき，c − uvZ 座標系を画像座標系と呼ぶ．このようなカメラモデルにおい
て，三次元空間上の点 M は，直線 OM とスクリーンとの交点である点 m に投
影される．
図 15: ピンホールカメラモデル
12
カメラパラメータのうち各カメラの焦点距離やスケール，歪みを表すパラメー
タを内部パラメータと呼び，各カメラの位置，向きを表すパラメータを外部パ
ラメータと呼ぶ．3 × 3 の内部パラメータ行列 A，3 × 4 の外部パラメータ行列
P を用いると，世界座標系での三次元座標 X = (x, y, z)T とカメラ座標系での
座標 x = (u, v)T との関係は式 (7) で表される．
(x)
(u)
λ
v
1
= AP
y
z
1
(7)
ここで，λ は実数である．
本研究ではまずカメラごとに内部キャリブレーションを行った後，推定した
内部パラメータを用いて全てのカメラの外部キャリブレーションを行う．以下，
本研究で用いる内部，外部それぞれのキャリブレーション手法について述べる．
図 16: キャリブレーションによく用いられるチェ
スパターン
3.2.2
内部キャリブレーション
先述の通り，内部パラメータは各カメラの焦点距離やズーム値，歪みといっ
た情報を持つ．内部パラメータ行列 A は以下の要素によって式 (8) のように構
成される．
( αu
A=
0
0
γ u0
α v v0
0 1
) ( f su
=
13
0
0
f k s u0
f s v v0
0
1
)
(8)
ここで，f は焦点距離，su , sv はそれぞれ u 方向，v 方向の画素サイズ，点 (u0 , v0 )
はスクリーンとカメラ光軸の交点の画像座標系での画素，ks は画像平面上の u
軸と v 軸の傾きを表すせん断係数である．
また，この内部パラメータ行列以外にもレンズによる歪みの影響を考慮する
必要がある．理想的なピンホールカメラにおいて点 N が点 nu = (uu , vu ) に投影
される場合でも，レンズによる歪みの影響で，点 nd = (ud , vd ) に投影されてし
まう．それぞれの方向のレンズ歪みを du , dv とすると，点 nu と点 nd の関係は
式 (9),(10) のようになる．
uu = ud + du
(9)
vu = vd + d v
(10)
さらに，du , dv は歪みパラメータ k1 , k2 を用いて式 (11),(12) のように計算される．
du = ud (k1 rd2 + k2 rd4 )
(11)
dv = vd (k1 rd2 + k2 rd4 )
(12)
ここで rd は式 (13) の通りである．
rd =
√
u2d + vd2
(13)
内部キャリブレーションは基本的にカメラセッティングを固定し，位置やサイ
ズが既知である参照物体を様々な方向から撮影しパラメータを獲得する手法が取
られている．今回はチェスボードを用いた Zhang のキャリブレーション手法 [10]
を用いる．この手法では，まず図 16 のような平面上に黒と白のチェスパターン
をプリントアウトしたものをカメラで観測し，その縦線と横線の交点 (コーナー)
を特徴点として検出する．この手法では，原理的に全てのコーナーが検出され
た画像が 3 枚以上あれば連立方程式を解くことで内部パラメータ αu , αv γ, u0 , v0
及び歪みパラメータ k1 , k2 を求めることができる．なお，コーナー検出の際，コ
ンピュータビジョン用のライブラリ，OpenCV[11] にすでに実装されているもの
を用いる．
3.2.3
外部キャリブレーション
先述の通り，外部パラメータは各カメラの位置関係や向きといった情報を持
つ．外部パラメータ行列 P はカメラ座標系から世界座標系への変換を表し，ま
14
ずカメラ座標系を回転させ，その姿勢を保ったまま平行移動させることで世界
座標系を得ることができる．P は回転行列 R と並進行列 T を用いて式 (14) の
ように表される．
P = [R|T]
(14)
ここで，座標系を x, y, z 軸周りにそれぞれ θx , θy , θz 回転させた場合の回転行
列 R は式 (15) のようになり，
R =
=
( cos θz
(
sin θz
0
r1 r2
r4 r5
r7 r8
−sinθz 0
cosθz 0
0
1
)
r3
r6
r9
)(
cos θy 0 sin θy
0
1
0
− sin θy 0 cos θy
)( 1
0
0
0 cos θx − sin θx
0 sin θx cos θx
)
(15)
x, y, z 軸方向にそれぞれ tx , ty , tz 平行移動させた場合の並進行列 T は式 (16) の
ようになる．
T=
( tx )
ty
tz
(16)
外部キャリブレーションについては Ueshiba らの手法 [12] を用いる．キャリブ
レーションしたい内部キャリブレーション済みのカメラ全てで同一の点を観測
し，各カメラでの画像内でのその点の座標の組を獲得することで外部パラメー
タ (外部パラメータ行列 P の要素 12 個) を求めることができる．今回は内部キャ
リブレーションで用いたチェスボード (図 16) を用いて同様にコーナー検出を行
い座標セットを多数獲得し，外部パラメータを求める．
3.3
三次元形状推定
本節では三次元形状の推定手法について述べる．赤色のレーザ光を用いてい
るので，チンダル現象下で入射光及び反射光の散乱光を観測すると図 17 のよう
に入射光および反射光の道筋が赤色の領域として得られる．全てのカメラでこ
のような観測画像を撮影し，それらの画像を用いて三次元空間での入射方向，反
射位置を求める．そのためにまず各画像上での入射方向，反射位置を推定する．
画像中でどの部分が入射方向，反射位置に相当するかの推定について，以下
に述べる性質を利用する．入射方向は入射光が一本の細いレーザ光であるため，
画像中で赤く細長い半直線として観測される．また，反射位置に関しては入射
15
光が到達する点であるため，入射光に相当する細長い直線上に存在する．
また，反射位置で反射した光は反射位置から遠ざかるにつれて減衰する．反
射位置を中心とした単位球と半径 r の球を考える．単位球の表面積 S の範囲を
通る光は半径 r の球において表面積 Sr2 の範囲を通るため，単位面積を通過す
る反射光の光量は反射位置からの距離 r の二乗に反比例する．よって，反射位
置から r 離れた位置での反射光の光量を実数 l を用いて
l
r2
と表すと，そのうち
カメラに向かう散乱光の割合が d の場合，カメラで観測される散乱光の光量は
dl
r2
となる．つまり，画像中では反射位置すなわち r が最小となる位置に相当す
る画素が周囲に比べて明るく観測され，反射位置から離れるにつれて反射光が
減衰することが言える．そのため，画像中の反射位置は入射光に相当する細長
い直線上で最も輝度値が高い点となる．以下，観測画像から以上で述べたよう
な入射方向及び反射位置を推定する手順を述べる．
まず画像から入射方向に相当する直線を推定することを考える．観測画像の
各画素は RGB それぞれ 8bit の値を持つ．しかし，入射光の領域を抽出するた
めには画像内でどの部分が赤くどの部分がそうでないかという情報さえあれば
よいので，その後の画像処理がしやすいように観測画像の赤成分に関する二値
化処理を施し，図 18 のような白と黒の二値画像を得る．この画像の白い領域
を，細長い入射光の領域とそれ以外の反射光の領域に分離したい．先ほど述べ
た通り，入射光の領域は細い半直線状となるため，画像処理の一つである収縮
処理を適当な回数かけると入射光の領域のみを消すことができる．これを利用
し，図 18 のような二値画像に収縮処理をかけた後，それと同じ回数だけ膨張処
理をかけると図 19 のように入射光の領域のみが消えた二値画像が得られる．図
18 に見られる入射光の領域が図 19 では消えているため，図 18,19 のような二値
画像の各画素 (u, v) の値をそれぞれ ia (u, v), ib (u, v) としたとき，式 (17) で表さ
れる ic (u, v) を各画素の値とした二値画像に直すと図 20 のようになり，入射光
の領域だけを残した二値画像を生成することができる．なお，式 (17) における
「·」は論理積，i は i の否定を表す．この二値画像に対して線分検出をかけるこ
とによりレーザ光に相当する直線を推定する．
ic (u, v) = ia (u, v) · ib (u, v)
(17)
さらに反射位置は今求めた入射光に相当する直線上で最も明るい点，すなわ
16
ち輝度値が最も大きい点とした．図 17 の画像から以上の情報を推定した結果を
図 21 に示す．
図 17: 観測画像
図 18: 二値化処理後の画像
図 19: 収縮膨張処理後の画像
図 20: 入射光領域を抽出した画像
図 21: 入射方向，反射位置推定結果
17
以上のようにして各カメラでの観測画像から推定した入射方向，反射位置を
用いて三次元空間での入射方向，反射位置を推定することができる．三次元空
間での反射位置 Ms = (x, y, z)T を推定する際の概念図を図 22 に示す．この図
にあるように，2 台のカメラ C1 , C2 の内部パラメータ行列をそれぞれ A1 , A2 と
し，外部パラメータ行列をそれぞれ P1 , P2 とし，画像内での反射位置をそれぞ
れ ms1 = (u1 , v1 )T , ms2 = (u2 , v2 )T とすると，式 (18),(19) のようになる．
λ1
λ2
( u1 )
v1
1
( u2 )
v2
1
(x)
= A 1 P1
y
z
1
(18)
(x)
= A 2 P2
y
z
1
(19)
この連立方程式を解くことで三次元空間での反射位置 Ms = (x, y, z)T 及び λ1 , λ2
を求めることができる．三次元空間での入射方向は，先程の反射位置 Ms と入
射光を表す半直線上の点 Me によって計算できる．つまり，Me の三次元座標を
求めることで三次元空間での入射方向を得ることができる．
図 22: 三次元空間での反射位置推定の概念図
Me の三次元座標を得るために 2 枚の観測画像での入射光上の対応点を求め
る．その求め方を以下に述べる．図 23 に示すように，まず，カメラ C1 の画像
上で入射光に相当する半直線上の反射位置と異なる点 me1 を選ぶ．その点のエ
ピポーラ線をカメラ C2 の画像に描き，エピポーラ線とカメラ C2 の画像上での
18
入射光に相当する半直線の交点 me2 を求める．この me1 と me2 が入射光上の対
応点となるので，この me1 , me2 に対して先程三次元空間での反射位置 Ms を求
めた際と同様にしてその点 Me の三次元座標を求める．こうして求めた Ms と
Me を持つベクトルを三次元空間での入射方向とする．Me 及び入射方向を推定
する際の概念図も図 24 に同様に示す．
図 23: 2 枚の観測画像での入射光上の対応点の求め方
図 24: 三次元空間での入射方向推定の概念図
19
第4章
4.1
チンダル現象下での BRDF 推定
チンダル現象下での BRDF 推定手法の概要
三次元形状を推定する際は画像中でどこが明るくてどこが暗いかといった輝
度値の大小の情報があれば十分であったが，BRDF を推定する際は輝度値の大
小だけでなく，入射光に対する反射光の光量の比を画像の輝度値から推定する
必要がある．画像で観測される輝度値は，カメラに直接入る反射光の他に，チ
ンダル現象による散乱光の影響を受けており，輝度値から直接 BRDF を推定す
るためには，BRDF と同時に散乱現象の影響も推定しなければならない．
BRDF に関しては 2 章で述べた通り，数式モデルを用いず，各反射方向への
反射率を求めることで推定する．一方散乱現象については，媒体中でどの方向
にどれほどの光が散乱し，どれほどの光が媒体に吸収され，どれほどの光がそ
のまま通過するかといった情報を数値で表した散乱パラメータを求めることで
推定する．
本研究ではその各反射方向への反射率及び散乱パラメータをパラメータ探索に
よって同時に推定する．具体的には，3.3 節で推定した入射光に関する情報と，与
えられた BRDF パラメータ及び散乱パラメータの下で各カメラで観測される画
像を光線追跡法を用いて生成する (以降，これをシミュレート画像と呼ぶ)．以下，
カメラ k での観測画像 Ikobs をシミュレートした画像を Iksim とし，その画素 (u, v)
での輝度値を観測画像の場合と同様に iksim (u, v) とする (u = 1, 2, · · · , Wsim ; v =
1, 2, · · · , Hsim ; シミュレート画像の解像度は Wsim × Hsim )．実際の観測画像と
それに対するシミュレート画像の誤差を計算する誤差評価関数を定義し，その
値が最小となるような各反射方向への反射率及び散乱パラメータを求める．
本章では以下，本研究における実物体の BRDF の推定手法について述べる．
まずチンダル現象によって生じる散乱現象に関して 4.2 節で述べ，シミュレー
ト画像の生成方法に関して 4.3 節で述べ，誤差評価関数に関して 4.4 節で述べ，
パラメータ探索アルゴリズムに関して 4.5 節で述べる．
4.2
散乱現象
4.1 節で述べた散乱パラメータは一般に 3 種類あると言われており，媒体の濃
度に依存するパラメータ 2 種類と媒体の種類に依存するパラメータ 1 種類に大
別される [13]．そのうち濃度に依存するパラメータには散乱係数と吸収係数が
20
ある．散乱係数とは散乱媒体中をある光量の光が単位距離を通過する際にその
光量のうち散乱する光量の割合であり，吸収係数とは媒体中をある光量の光が
単位距離を通過する際に媒体中で吸収される光量の割合である．散乱係数 β ，吸
収係数 µ とおくとこれらは濃度により．0 ≤ β ≤ 1; 0 ≤ µ ≤ 1 の範囲で変化す
る．図 25 のように光量 L の光が単位距離を通過するとき，散乱する光の光量
Ls ，微粒子に吸収される光量 La は式 (20),(21) のように表される．
Ls = βL
(20)
La = µL
(21)
結果，散乱も吸収も起こらず直進する光量 Lt は式 (22) のように表される．
Lt = (1 − β − µ)L
(22)
図 25: 散乱媒体による散乱現象
散乱光 Ls の散乱方向に関する分布は，入射光 L の進行方向と散乱方向のな
す角 θ と媒体の種類に依存し，パラメータ g によって式 (23) の形で与えられる．
P (g, θ) =
1
1 − g2
4π (1 + g 2 − 2gcosθ)3/2
21
(23)
これを用いると，入射光との角度が θ となるある方向に散乱する光の光量 Ls (θ)
は式 (24) のように表される．
g の値によって散乱の振る舞いは変化する．g = 0 のときは全ての方向に同じ
光量の散乱光が発生する等方散乱となり，g > 0 のときは元々の光の進行方向
と同じ方向付近で散乱光の光量が多くなる前方散乱となり，g < 0 のときは元々
の光の進行方向と逆の方向付近で散乱光の光量が多くなる後方散乱となる．図
26 に示すように，g の値が大きくなると進行方向付近の散乱光の割合が大きく
なる．
Ls (θ) = βP (g, θ)L
(24)
図 26: 散乱関数パラメータ g の散乱光分布への影響
Narasimhan らによるとほとんどの媒体は前方散乱に分類される [13]．その例
として，媒体として煙を用いたものを図 27 に示す．図 27 は，計測装置内でレー
ザを照射した際に，レーザの進行方向に対してその散乱光を前方から観測した
場合 (左図) と後方から観測した場合 (右図) の画像を示している．両画像を比較
すると光の前方から観測した方が光量の多い散乱光を観測することができ，レー
ザの進行方向とカメラの位置関係によって観測される散乱光の強さが異なるこ
とが分かる．従って，画像上で観測される散乱光の輝度値から入射光および反
射光の光量を推定するためには，本節の冒頭でも述べた通り，輝度値だけでな
く散乱現象による影響を考慮する必要がある．また，散乱現象による影響を考
慮するためには入射光および反射光の進行方向とカメラの位置関係も加えて考
慮する必要がある．
22
図 27: チンダル現象下でレーザ光の前方及び後方から散乱光を
観測した様子
4.3
4.3.1
シミュレート画像生成
光線追跡
本節では，光線追跡法を用いたシミュレート画像生成方法について述べる．ま
ず光線追跡の概要を述べる．レーザ照射装置から照射された入射光及び物体表
面で反射した後の反射光は媒体中で散乱し，各カメラに向かう．その散乱光の
光量を計算する．
散乱媒体中で散乱した光が各カメラに到達する．ここで，カメラのスクリー
ン上のある画素 m に一次散乱光が到達する様子を図 28 に示す．この図におい
て画素 m に到達する散乱光の光量は，カメラ中心と画素 m とを結ぶ視線上の
点で散乱してカメラに届く光の光量を全て足しあわせたものとなる．
シミュレート画像の各画素に到達する光の光量を計算するためには，三次元
空間の各点における散乱光の光量を計算する必要がある．そのために三次元空
間を単位空間に分割し，入射光および反射光が各単位空間を通過する際に吸収
される光の光量，散乱される光の光量，通過する光の光量をそれぞれ計算する．
4.3.2
単位空間の定義
本節では，三次元空間の単位空間分割をどのように行うかについて述べる．通
常よく用いられる立方格子状の空間分割を行った場合，一般には一つの単位空
間を通過する際に光が進む距離は一定とならない．例えば図 29 に示した立方格
23
図 28: ある画素に届く散乱光
子状の空間分割の断面図において，入射光が水色の単位空間を通過する際に光
が進む距離は黄緑色の単位空間を通過する際に光が進む距離に比べて短くなる．
散乱係数は単位距離を通過する際に散乱される光の割合であるため，このよう
に単位空間ごとに光の進む距離が変動する場合は散乱光の光量などの計算が複
雑になってしまう．
そのため，本研究では反射位置を原点とした球座標系の空間分割を用いる．入
射光及び反射光は原点である反射位置を通るため，図 30 に示した球座標系の空
間分割の断面図において，入射光が水色の単位空間を通過する際も黄緑色の単
位空間を通過する際も光が進む距離が同じになるように，どの単位空間を通過
する際も光が進む距離が一定となり計算が簡単になる．この球座標系の空間分
割の詳しい定義を以下で述べる．
図 29: 立方格子状の空間分割
24
図 30: 球座標系の空間分割
まず反射位置を中心とした直交座標系をとる．反射位置 X0 = (x0 , y0 , z0 )T を
原点，z 軸を反射位置を含む接平面の法線方向，x 軸は入射方向を xy 平面 (z 軸
と直交し原点を通る平面) に投影した方向とし，右手直交座標系になるように y
軸をとる．この座標系を反射位置中心直交座標系と呼ぶ．このとき，反射位置中
心直交座標系の座標 X′ = (x′ , y ′ , z ′ )T を世界直交座標系の座標 X = (x, y, z)T に
変換する際の回転行列を R′ , 並進ベクトルを T′ とするとこれらの関係は式 (25)
のようになる．
X = R′ X′ + T′
(25)
さらにこの反射位置中心直交座標系を反射位置中心球座標系 (r, θ, ϕ) に変換
する．その変換は式 (26),(27),(28) で表される．以下，r, θ, ϕ で表される球座標
系単位空間を V (r, θ, ϕ) と表記する．
r=
√
x′2 + y ′2 + z ′2
z′
r
y′
ϕ = arctan ′
x
θ = arcsin
4.3.3
(26)
(27)
(28)
シミュレート画像生成アルゴリズム
シミュレート画像は 4.1 節で述べた通り，与えられた各反射方向への反射率及
び散乱パラメータや 3.3 節で推定した三次元空間での入射方向と反射位置から
25
光線追跡により生成する．光線追跡によるシミュレートの対象となるのは，反
射位置に向かって進む入射光及び反射位置から様々な方向に向かって進む反射
光である．3.1 節で述べた通り，多重散乱は考慮しない．以下にシミュレート画
像生成アルゴリズムを述べる．このアルゴリズムにおいて，a → b は a を b に
更新することを表す．
アルゴリズム
このアルゴリズムでは，入射光が反射位置に向かって単位空間を
一つずつ進み，反射位置で反射し，反射位置から様々な方向に向かって単位空間
を一つずつ進む間に各単位空間で散乱してカメラに届く光の光量を計算し，各
カメラの各画素に到達する光量から各カメラでのシミュレート画像を生成する．
1. 光線追跡によるシミュレートを開始する入射位置を決定する．球座標系の
空間分割を反射位置からの距離 r が rmax 以内の単位空間で構成される空間
としたとき，入射方向が θi , ϕi ならば，V (rmax , θi , ϕi ) が入射位置の単位空
間となる．また，この時点での光の光量 Lnow を L，単位空間を Vnow とし，
シミュレート画像の各画素の輝度値 ikobs (u, v) はすべて 0 に設定する．
2. 光が単位距離進む際の散乱現象を計算するために，反射位置に向かって単
位距離 (単位空間一つ分) 進む．Vnow が V (r, θi , ϕi ) であるとするとこれを
V (r − 1, θi , ϕi ) に更新する．
3. Vnow に入射した光量 Lnow のうち，各カメラ k へ散乱する光量 Lks を散乱係
数 β ，散乱関数パラメータ g を用いて計算する．元々の光の進行方向と散
乱位置からカメラ k に向かう方向の間の角度を θk とすると Lks は以下の式
(29) と式 (24) により計算できる．各カメラに対してこれを計算する．
Lks = Ls (θk )
(29)
4. 現在の反射位置中心球座標を世界直交座標に変換し，カメラキャリブレー
ションの結果から式 (7) を用いてその位置が各カメラのどの画素 (u, v) に投
影されるかを計算する．式 (30) の通りにその画素に到達した光量 Lks を加
える．
iksim (u, v) → iksim (u, v) + Lks
(30)
5. 現在の単位空間を通過する光量を散乱係数 β ，吸収係数 µ を用いて式 (22) の
ように計算する．これが次の単位空間に到達する光量となるので，Lnow →
Lt と更新を行う．
26
6. 2-5 を反射位置に相当する単位空間に到達する (r = 0) まで繰り返す．反射
位置に到達した光の光量を L′ とする．
7. 各反射方向 ((90/∆θ) × (360/∆ϕ) 通り,∆θ, ∆ϕ は反射方向の刻み幅) ごとに
8-11 を行う．
8. 与えられた反射率によりその方向への反射光の初期光量を計算する．このと
きの反射方向 θr = i, ϕr = j(i = 0◦ , ∆θ◦ , · · · , (90−∆θ)◦ , j = 0◦ , ∆ϕ◦ , · · · , (360−
∆ϕ)◦ ) に対する反射率 fr (i, j) を用いて Lnow → fr (i, j)L′ と更新を行う
(0 ≤ fr (i, j) ≤ 1)．
9. 2 と同様に今度は反射位置から単位距離 (単位空間一つ分) 離れる方向に進
む．Vnow が V (r, θr , ϕr ) であるとするとこれを V (r + 1, θr , ϕr ) に更新する．
10. 3-5 と同様の計算を行う．
11. 9-10 を反射光が最も外側の単位空間 V (rmax , θr , ϕr ) に到達するまで繰り返す．
12. 各画素の輝度値 iksim (u, v) から各カメラでのシミュレート画像 Iksim を生成
する．
なお，3. で計算した各単位空間で散乱しカメラ k に届く光量 Lks を 4. のように
足しあわせて各画素の輝度値 iksim (u, v) を計算し，13. でシミュレート画像 Iksim を
生成する際，坂本らが開発したオープンソースライブラリ，Kyoto Visualization
System(KVS)[14] を用いた．
4.4
評価関数の定義
4.3.2 節のアルゴリズムで生成されたシミュレート画像と観測画像を誤差評価
関数を用いて比較し，それを最小にする各反射方向への反射率及び散乱パラメー
タを探索する．本節では本研究での誤差評価関数の詳細を述べる．
観測画像 Ikobs とシミュレート画像 Iksim の解像度が同じ (Wobs = Wsim , Hobs =
Hsim ) 場合，両者の誤差は各画素同士の輝度値の差を足しあわせた式 (31) で表
される E k とする．
k
W
H
∑
∑
E =
u=1
|iksim (u, v) − ikobs (u, v)|
(31)
v=1
しかし，観測画像におけるシャッタースピードや入射光の光量によって光の絶
対的な明るさが変化することを考慮し，シミュレート画像の各画素の輝度値を
正規化する．具体的には，正規化によって観測画像とシミュレート画像の各画
27
素の輝度値の和が同じになるようにする．観測画像において強い光が観測され
る部分は最大値 (今回は 255) に達していることがあるため，単純にシミュレー
ト画像の各画素を定数倍することはできない．そこで，観測画像の輝度値 255
に対応するシミュレート画像の輝度値 ikstd を決定し，シミュレート画像におい
てそれ以上の輝度値を持つ画素は 255 に丸めるという正規化補正を行う．
その処理手順を述べる．まず，観測画像の各画素の輝度値の和 Sumkobs を式
(32) のように計算する．
Sumkobs =
W
H
∑
∑
u=1
次に，ikstd
ikobs (u, v)
(32)
v=1
をもとにシミュレート画像の各画素の輝度値を式 (33) に従って補
正する．
i′k
sim (u, v) = min(
255 · iksim (u, v)
ikstd
, 255)
(33)
こうして補正したシミュレート画像の各画素の輝度値の和 Sum′k
sim を式 (34)
k
のように計算し，式 (35) で計算される Sumkobs と Sum′k
sim の差 Sub が最小とな
る ikstd を求める．
Sum′k
sim =
W
H
∑
∑
u=1
i′k
sim (u, v)
(34)
v=1
k
Subk = |Sum′k
sim − Sumobs |
(35)
求めた ikstd を用いてシミュレート画像のすべての輝度値を補正し，その画像
′k
k
I′k
sim を得る．Isim と Iobs の各画素の輝度値の差を足しあわせることで，式 (36)
のように E k を求める．各カメラに置いて E k を計算し，E k を式 (37) のように
全てのカメラに関して足しあわせた値 E を評価関数の値とする．
Ek =
W
H
∑
∑
u=1
k
|i′k
sim (u, v) − iobs (u, v)|
v=1
28
(36)
K
∑
E=
Ek
(37)
k=1
4.5
パラメータ探索による BRDF 推定
本節ではパラメータ探索の方法について述べる．式 (35) で示した誤差評価関数
E を最小とする各反射方向への反射率及び散乱パラメータを Downhill Simplex
Method[15] を用いて獲得することにより BRDF を推定する．
このアルゴリズムを用いて探索を行う際に，適切な初期値を与えることによっ
て探索にかかる時間を短縮し，誤った局所解に陥る可能性を減らすことができ
る．最終的に探索を行うパラメータは (90/∆θ) × (360/∆ϕ)+3 種類で膨大とな
るので，最初はパラメータの数を少なくした状態で探索を行い，その結果を初
期値とし改めて探索を行う．パラメータの数を少なくするために，式 (2) で示
される Lambert モデルと Torrance-Sparrow モデルを利用する．
このモデルの 4 つのパラメータ拡散反射係数 kd , 鏡面反射係数 ks , 表面粗さ係
数 σ, 屈折率 n によって式 (38) の BRDF を式 (39) のように表現することができ
る．ここで，式 (38) の θi , ϕi が式 (39) の L に，式 (38) の θr , ϕr が式 (39) の V に
それぞれ対応する．なお，式 (39) の N は法線方向なので固定である．以降，こ
のモデルのパラメータを推定するという前提で手法の説明を行うが，他の数式
モデルを用いてもパラメータ数が変更するだけで基本的には同じである．
fr = f (θi , ϕi , θr , ϕr )
fr = kd (N · L) +
ks
exp
N·V
(
−
)
2
α
2σ 2
(38)
·G·F
(39)
この反射モデルを用いて，本研究では大きく 2 つのフェーズに分けてパラ
メータ空間の探索を行う．1 つ目は各反射方向への反射率を Lambert モデル及
び Torrance-Sparrow モデルで近似し，パラメータを Torrace-Sparrow の反射モ
デルのパラメータ 4 種類と散乱パラメータ 3 種類の計 7 種類として探索を行い，
各反射方向への反射率及び散乱パラメータの適切な初期値を求めるフェーズで
ある．数種類のパラメータで近似した反射モデルを用いることからこれをパラ
メトリックフェーズと呼ぶ．
2 つ目はパラメトリックフェーズの結果のうちの Lambert モデル及び Torrance29
Sparrow モデルのパラメータ 4 種類から計算される (90/∆θ) × (360/∆ϕ) 方向の
反射率 fr (i, j) を初期値とし，パラメトリックフェーズの結果のうちの散乱パラ
メータ 3 種類をそのまま加えた (90/∆θ) × (360/∆ϕ)+3 種類のパラメータ空間
の探索を行うフェーズである．これを先程のパラメトリックフェーズに対して
ノンパラメトリックフェーズと呼ぶ．このように 2 つのフェーズに分けること
で適切な初期値を与えることができるため，探索にかかる時間を短縮し，誤っ
た局所解に陥る可能性を減らすことができる．
30
第5章
実験
提案手法の有効性を検証するための実験を行った．まず，本実験で使用する
計測装置について 5.1 節で述べる．次に，誤差評価関数の振る舞いを測定し，最
適解への収束性と解の一意性を評価した．これを 5.2 節で述べる．最後に，実
際の観測画像に対して提案手法によるパラメータ探索を行い，三次元形状，各
反射方向への反射率，散乱パラメータの推定を行った．これを 5.3 節で述べる．
5.1
計測装置
3.1 節で述べた観測環境として以下に述べる計測装置を用いた．今回の実験で
は，図 31 に示すようなアクリル板で囲まれた 80cm 立方の密閉された空間を用
い，その中にカメラ 5 台を固定し，密閉空間内の散乱媒体の濃度を均一に保つ
ためのファンを設置し，底面の中央に計測対象を置いた．密閉空間外では，プ
ログラムにより制御可能な赤色レーザポインタを取り付けたレーザ照射装置 (図
32) 及び専用の液体を気化させた煙を散乱媒体として密閉空間に充満させるた
めのスモークマシン (図 33) がある．スモークマシンから発生させた煙を密閉空
間内に充満させチンダル現象を生じさせた．
今回用いるカメラに関しては全てデジタル一眼レフカメラ「Nikon D40」を
用いており，3008 × 2000 の高解像度となっている．あらかじめ，ステレオ計測
が可能となるように様々な方向から観測できるようにカメラを配置し，底面の
中央に置いた計測対象が画像の中心に映るようにカメラの向きやズームを調節
し，カメラキャリブレーションを行った．
図 31: 計測装置
図 32: レーザ照射装置
31
図 33: スモークマシン
5.2
誤差評価関数の正当性の評価
本節では，誤差評価関数の振る舞いを測定し，最適解への収束性と解の一意
性を評価した．
具体的な実験方法としては，まずあるパラメータ (以降これを正解パラメー
タと呼ぶ) からシミュレート画像を生成し，これを入力画像とした．そして，そ
のパラメータを変化させたシミュレート画像との誤差を誤差評価関数によって
求めた．パラメータ空間が広くその全体について調べることは難しいため，正
解パラメータ付近でパラメータを変化させながらシミュレート画像を生成し正
解パラメータから生成したシミュレート画像との誤差を求め，誤差評価関数が
正解パラメータにおいて極小となることを示した．
なお，正解パラメータは，散乱係数 β = 1.0 × 10−4 ，吸収係数 µ = 1.0 × 10−4 ，
散乱関数パラメータ g = 5.0 × 10−1 ，拡散反射係数 kd = 3.0 × 10−5 ，鏡面反射
係数 ks = 5.0 × 10−3 ，表面粗さ係数 σ = 1.0 × 10−1 ，屈折率 n = 1.5 とした．こ
のパラメータにより生成したカメラ 5 台分のシミュレート画像を図 34 に示す．
5.2.1
パラメータ 1 種類を変化させた場合
本節では，正解パラメータ 7 種類のうち 6 種類を正解パラメータに固定し，残
りの 1 種類を変化させ，そのパラメータから生成したシミュレート画像と正解
パラメータから生成したシミュレート画像での誤差評価関数の変化を調べた．
図 35-41 にそれぞれ散乱係数，吸収係数，散乱関数パラメータ，拡散反射係数，
鏡面反射係数，表面粗さ係数，屈折率を正解パラメータから変化させたときの
誤差評価関数の変化を示す．各パラメータを変化させる範囲はそれぞれのパラ
メータの定義域であり，その定義域は表 1 のようになっている1) ．なお，表面粗
さ係数と屈折率の上限に関してはそれぞれ馬場らの研究 [16] と屈折率，フレネ
ル値一覧表 [17] を参考に決定した．
1)
β と µ には β + µ ≤ 1.0 という制約もある．
32
図 34: 正解パラメータでのシミュレート画像
上限
下限
β
µ
g
kd ks
σ
1.0 1.0 1.0 1.0 1.0 1.0
0.0 0.0 -1.0 0.0 0.0 0.0
表 1: 各パラメータの定義域
n
5.0
1.0
図 35: 散乱係数のみの変化
パラメータの軸の赤い丸で囲った値がそれぞれの正解パラメータであり，こ
れらの図に示すように，7 種類のパラメータのうちどれか 1 種類を変化させた
場合，正解パラメータにおいて誤差評価関数の値が極小となることが分かる．
5.2.2
パラメータ 2 種類を変化させた場合
次に，7 種類のパラメータのうち 2 種類を変化させ，そのパラメータから生成
したシミュレート画像と正解パラメータから生成したシミュレート画像での誤
33
図 36: 吸収係数のみの変化
図 37: 散乱関数パラメータのみの変化
図 38: 拡散反射係数のみの変化
差評価関数の変化を調べた．今回誤差評価関数の変化を調べる範囲は，2 種類の
パラメータ par1 , par2 を 5.2.1 節においてそれぞれ単独で変化させた場合の誤差
評価関数が 1.0 × 106 以内となる範囲がそれぞれ min1 ≤ par1 ≤ max1 , min2 ≤
par2 ≤ max2 であるとすると，min1 ≤ par1 ≤ max1 かつ min2 ≤ par2 ≤ max2
を満たす範囲とする．7 種類のパラメータから 2 種類を選ぶ組み合わせは 7 C2 通
34
図 39: 鏡面反射係数のみの変化
図 40: 表面粗さ係数のみの変化
図 41: 屈折率のみの変化
りあり，そのうち表面粗さ係数と屈折率を変化させた様子を図 42 に，散乱係数，
吸収係数を変化させた様子を図 43 に，散乱係数，拡散反射係数を変化させた様
子を図 44 にそれぞれ示す．
このグラフでは色によって誤差評価関数の値の大きさを表しており，青→緑
→赤と変化するに連れて誤差評価関数の値が大きくなっていく．両パラメータ
の軸の赤い丸で囲った値がそれぞれの正解パラメータであり，これらの図に示
35
すように，2 種類のパラメータを変化させた場合でも正解パラメータにおいて
誤差評価関数の値が極小となることが分かる．この誤差評価関数を極小にする
ことで各反射方向への反射率および散乱パラメータを推定することができる．
よってこの誤差評価関数の定義が正当なものであると言える．
図 42: 表面粗さ係数と屈折率の変化
図 43: 散乱係数と吸収係数の変化
36
図 44: 散乱係数と拡散反射係数の変化
5.3
観測画像を用いた三次元形状，BRDF パラメータ，散乱パ
ラメータ推定
本節では，実際の観測画像に対して提案手法によるパラメータ探索を行い，三
次元形状，各反射方向への反射率，散乱パラメータの推定を行う．今回用いる
対象物体として，従来計測が困難であった黒色物体及び金属物体を選択した．
5.3.1
黒色物体
黒色物体として用いるのは図 45 に示す黒い手帳である．この物体にチンダ
ル現象下で赤色レーザポインタを照射した様子を 5 台のカメラで観測した画像
を図 46 に示す．これらの画像に対して提案手法による入射光，反射位置推定を
行う．
図 45: 黒い手帳
37
図 46: 観測画像
図 46 の各画像から入射方向及び反射位置を推定した結果を図 47 に示す．こ
の図でも金属物体の実験と同様に入射方向を緑の線分で，反射位置を緑の円で
示している．この図から，提案手法によりそれぞれ正確に推定できていること
が分かる．その後三次元空間での入射方向及び反射位置をステレオ視により計
算した．なお，図 46 の画像を見ると，黒色物体の場合は反射光がほとんど発生
していないことが分かる．そのため，この物体に関しては各反射方向への反射
率および散乱パラメータのパラメータ探索は割愛する．
5.3.2
金属物体
金属物体として用いるのは図 48 に示す金属球である．この物体にチンダル現
象下で赤色レーザポインタを照射した様子を 5 台のカメラで観測した画像を図
49 に示す．これらの画像に対して提案手法による入射光，反射位置推定及びパ
ラメータ探索を行った．
三次元形状推定
図 49 の各画像から入射方向及び反射位置を推定した結果を図
50 に示す．この図では入射方向を緑の線分で，反射位置を緑の円で示す．この
図から，提案手法によりそれぞれ正確に推定できていることが分かる．その後
三次元空間での入射方向及び反射位置をステレオ視により計算した．
合計 296 フレーム撮影した結果から上記のようにして 253 点の反射位置を獲
得した．正解の形状に獲得した三次元位置の集合をプロットし，二方向から見
た図が図 51(a),(b) である．図の中で青色の点が獲得した三次元位置で，赤色の
点の集合が正解の形状である．
38
図 47: 入射方向及び反射位置推定結果
獲得した点がどれだけ正解の金属球の表面に近いかを示す尺度として，金属
球の中心からの距離を用いる．金属球の半径を r，推定した三次元位置の金属
球の中心からの距離を r′ とし，それぞれの三次元位置について |r − r′ | を計算し
て評価する．その度数分布を表 2 に，ヒストグラムを図 52 に示す．
図 51 にあるように，提案手法により獲得した点群により大まかな球体の形状
を獲得できた．また，表 2，図 52 にあるように，獲得した三次元位置と正解の
39
図 48: 金属球
図 49: 観測画像
球の表面との距離 (誤差) は概ね 3mm 以内に収まっている．また誤差の平均は
1.078mm であった．
各反射方向への反射率及び散乱パラメータ推定
図 50 の結果を元に観測画像か
ら各反射方向への反射率及び散乱パラメータをパラメータ探索により推定した．
ただし今回は実装の都合上観測画像 Ikobs より解像度が小さいシミュレート画像
Iksim しか生成することができなかったため，誤差評価関数を求める前にまず観
測画像を縮小しシミュレート画像と解像度を合わせる処理を行った．
ここで Ikobs を縮小した後の観測画像を I′k
obs とし，その画像中の画素 (u, v) での
k
′k
輝度値を i′k
obs (u, v) とする．今回は Iobs の画素 (u, v) に対応する画像 Iobs 中の領域
内の全画素の平均輝度値を i′k
obs (u, v) とする．例えば，図 53 中の左の大きい画像
中の緑色の部分は右の小さい画像中の赤色の画素に対応する領域である．観測
40
図 50: 入射方向及び反射位置推定結果
画像 Ikobs とシミュレート画像 Iksim の解像度の比 RW , RH をそれぞれ式 (40),(41)
で表し，
Wobs
Wsim
Hobs
RH =
Hsim
RW =
(40)
(41)
′k
I′k
obs の輝度値 iobs (u, v) は式 (42) のように対応する画素の輝度値の平均をとる．
41
(a)
(b)
図 51: 金属球における形状の推定結果
表 2: 金属球における反射位置の誤差 |r − r′ | の度数分布
誤差
頻度 累積% 平均誤差
0-1mm
1-2mm
2-3mm
3-4mm
4-5mm
5-6mm
6-7mm
7-8mm
8-9mm
9-10mm
10mm 以上
175
45
18
5
3
3
1
0
0
0
2
69.44%
87.30%
94.44%
96.43%
97.62%
98.81%
99.21%
99.21%
99.21%
99.21%
100.00%
1.078mm
全ての画素に対してこの計算を行い，縮小後の観測画像 I′k
obs を生成する．
∑
RW (u+1)
1
i′k
obs (u, v) =
RW RH
∑
RH (v+1)
ikobs (u′ , v ′ )
(42)
u′ =RW u+1 v ′ =RH v+1
′k
その後 I′k
obs と Isim の誤差を先述の通りに計算し，E を求める．
以下，パラメトリックフェーズ，ノンパラメトリックフェーズの順にパラメー
タの探索結果を述べる．
パラメトリックフェーズ
まず 7 種類のパラメータの初期値を決める．今回は
誤差評価関数の値が小さい初期値を使うことで探索時間を削減することを考え，
予め 7 種類のパラメータをランダムに変化させながらそれぞれの誤差評価関数の
値を計算したデータを多数取得し，その中で誤差評価関数の値が低いパラメー
タをいくつか選びそれらを初期値とした．これらの初期値からパラメータ探索
を開始した．その結果を表 3，4 に示す．表 3 は各パラメータの初期値を pars
42
図 52: 金属球における反射位置の誤差 |r − r′ | のヒスト
グラム
図 53: 縮小前と縮小後の画素の対応関係
とした場合の探索終了時のパラメータ pare をまとめたものである．便宜上，そ
れぞれの初期値に探索 ID を割り振る．また，表 4 は各探索 ID の探索における
pars での誤差評価関数の値 Es ，pare での誤差評価関数の値 Ee をまとめたもの
である．
表 3 に注目すると，様々な初期値から探索を行った結果，必ずしも同じパラ
メータに収束するわけではなく，初期値に依存した局所解にそれぞれ収束して
いることがわかる．その中でも，散乱関数パラメータ g ，表面粗さ係数 σ ，屈折
率 n に関しては各探索の結果において分散が小さくなっており，おおよそ正し
い値に収束していると考えられる．一方，散乱係数 β ，吸収係数 µ，拡散反射係
43
表 3: 探索開始時と終了時の各パラメータの値
探索 ID
1
2
3
4
pars
pare
pars
pare
pars
pare
pars
pare
β(10−3 )
0.42
2.59
0.99
19.44
0.66
15.16
0.36
0.00
µ(10−3 ) g(10−1 )
8.25
6.90
18.85
3.52
8.77
6.97
1.13
3.77
8.76
7.24
0.00
5.75
8.47
7.13
20.46
4.03
kd (10−3 )
0.57
3.11
2.19
4.25
0.88
2.27
0.65
4.10
ks (10−1 )
2.27
6.79
3.20
9.22
2.19
4.09
7.20
9.18
σ(10−2 )
4.45
2.20
2.60
2.24
5.87
2.50
2.48
2.21
n
2.20
4.99
1.84
5.00
2.29
5.00
2.13
4.90
表 4: 誤差評価関数の値
探索 ID
1
2
3
4
Es (106 )
2.07748
2.81397
2.18384
2.03780
Ee (106 )
1.97057
1.96395
1.98469
1.95684
数 kd ，鏡面反射係数 ks に関してはその分散が大きくなっている．ここで，探索
ID1 の pare でのシミュレート画像及び探索 ID2 の pare でのシミュレート画像を
それぞれ図 54，図 55 に示す．これらの画像について，いずれも観測画像に見ら
れる強く鋭い鏡面反射をほぼ再現できた画像となっている．
また，各探索 ID において推定した Lambert モデル及び Torrance-Sparrow モ
デルのパラメータを用いて環境マッピングによるレンダリングを行った．その
結果を図 56 に示す．レンダリング結果としては，計測対象である金属球におい
て推定したパラメータにおける鏡面反射係数 ks が拡散反射係数 kd に比べてか
なり大きいことや，表面粗さ係数が小さいという特徴を反映した，滑らかで金
属光沢のある表面反射特性を再現することができた．
しかし，ID1 と ID2 のシミュレート画像を比較すると類似した画像となって
おり，表 4 に示す通り誤差評価関数の値にも大きな差は見られない．
このように異なるパラメータで同様のシミュレート画像が生成される原因と
しては 2 つ考えられる．一つは入射光の画像上での輝度値が輝度値の上限 (255)
に達してしまっていることが挙げられる．その様子は図 57 に示す水色で囲んだ
部分などで確認することができる．この場合，拡散反射係数 kd と鏡面反射係数
ks との比は得ることができるが，kd と ks の絶対的な値を得ることはできない．
これが，kd と ks の推定に影響を与えていると考えられる．これを解決するため
には，入射光の画像上での輝度値が上限に達さないよう画像のダイナミックレ
44
図 54: 探索 ID1 の pare でのシミュレート画像
ンジを広くし，現状の 8bit 画像ではなく 16bit 画像などのさらに階調の多い画
像を用いるなどの工夫が必要となる．
もうひとつの原因としては，複数のパラメータを適度に変化させると物体表
面での反射や散乱媒体中の散乱現象単体としては変化が生じるが，それぞれの
変化の相殺によりほぼ同じシミュレート画像が生成されるということが考えら
れる．具体的には，以下の式 (43),(44) の条件を満たすパラメータの二種類組み
合わせ (β0 , µ0 , g, kd0 , ks0 , σ, n 及び β1 , µ1 , g, kd1 , ks1 , σ, n) から得られる二枚のシ
ミュレート画像はとほぼ同じ画像となる．
45
図 55: 探索 ID2 の pare でのシミュレート画像
β0 + µ0 = β1 + µ1
(43)
kd
ks
β1
= 0 = 0
β0
kd 1
ks1
(44)
そこで，この 7 種類パラメータのうち，散乱係数 β を 7.0 × 10−3 , 吸収係数
µ を 5.0×10−3 , 散乱関数パラメータ g を 6.5×10−1 に固定した上で残りの 4 種類
のパラメータの探索を行い，誤差評価関数の収束の一意性が高まるかどうかを
検証した．その結果を表 3，4 と同様の形式で表 5，6 に示す．表 3，4 に比べる
と，拡散反射係数 kd の分散が小さくなっていることが分かる．しかし，依然と
46
図 56: 環境マッピングによるレンダリング結果
して ks の分散が大きいままとなり，推定値が初期値の選び方に左右されている
ことが分かる．鏡面反射係数 ks に関しては今回のような計測対象の場合拡散反
射係数 kd が鏡面反射係数 ks に対して無視できるほど小さいためであると考え
られる．
初期値によって得られるパラメータが異なるもう一つの原因として，元々の
観測画像にはノイズが含まれていることが多く，正しいパラメータ以外にも極
小値が多数存在し，そのため，適切な初期値を選択しなければ正しいパラメータ
に収束せず，結果として誤った局所解に陥ってしまうということが考えられる．
今回の誤差評価関数では，観測画像とシミュレート画像の各画素の輝度値の
47
図 57: 観測画像中で輝度値の上限に達した明るい部分の例
表 5: 探索開始時と終了時の各パラメータの値
探索 ID
1
2
3
4
pars
pare
pars
pare
pars
pare
pars
pare
β(10−3 )
7.00
7.00
7.00
7.00
7.00
7.00
7.00
7.00
µ(10−3 ) g(10−1 )
5.00
6.50
5.00
6.50
5.00
6.50
5.00
6.50
5.00
6.50
5.00
6.50
5.00
6.50
5.00
6.50
kd (10−3 )
0.57
0.52
2.19
0.50
0.88
0.75
0.65
1.03
ks (10−1 )
2.27
2.13
3.20
4.50
2.19
2.88
7.20
6.69
σ(10−2 )
4.45
2.55
2.60
2.53
5.87
2.55
2.48
2.55
表 6: 誤差評価関数の値
探索 ID
1
2
3
4
Es (106 )
2.07318
2.71582
2.19506
2.00392
Ee (106 )
1.99099
1.99252
1.99171
1.99138
差を単純に足しあわせたが，誤差を計算する画素の輝度値に代わって周囲の画
素との重みつき和を用いるなど，この誤差評価関数を観測画像中のノイズに対
して頑健となるように設計すればこの局所解に陥りやすい問題に対処できると
考えられる．また，パラメータのうち β, µ, g を固定すると kd の推定結果の分散
が小さくなったように，7 つのパラメータを同時に探索するのではなく，その
一部を固定した上で探索を行えば，より誤差評価関数の収束の一意性が高まる
48
n
2.20
3.05
1.84
2.04
2.29
3.17
2.13
2.33
と考えられる．例えば，BRDF パラメータが既知の物体を利用してまず散乱パ
ラメータを推定し，そのパラメータを用いて BRDF パラメータが未知の物体の
BRDF パラメータを推定するという方法が挙げられる．また，kd , ks の比だけ
でなく絶対的な値を得るためには，ダイナミックレンジが広く階調の多い画像
を用いることが必要となる．
ノンパラメトリックフェーズ
次にパラメトリックフェーズの結果を用いて数
式モデルを用いない BRDF を各反射方向への反射率として推定した．パラメト
リックフェーズで様々なパラメータに収束したもののうち，表 3 の探索 ID1 の
結果を用いた．反射方向の刻み幅 ∆θ, ∆ϕ はともに 3 とした．
表 7 にノンパラメトリックフェーズの初期値 (=パラメトリックフェーズの推
定値)pars 及び推定値 pare における誤差評価関数の値 (それぞれ Es , Ee とする)
を挙げる．また，初期値 pars , pare でのシミュレート画像を図 58，59 に示す．こ
れらの誤差評価関数の値や画像に関して，初期値 pars と推定値 pare を比較し
た．特に誤差が改善されたカメラ 5 の観測画像，初期値 pars でのシミュレート
画像，推定値 pare でのシミュレート画像の全体図および拡大図を図 60 に示す．
拡大図に示された黄緑の楕円の内部を見ると，観測画像では鏡面反射光以外に
若干強い反射が見られる．パラメトリックフェーズではこれらの反射を再現で
きていないのに対し，推定値 pare での結果では，これを完全ではないものの再
現できている．
表 7: 誤差評価関数の値
探索 ID
1
Es (106 )
1.97057
Ee (106 )
1.93412
ノンパラメトリックフェーズにおける誤差の改善を示したのが図 61 である．
この図では，観測画像，パラメトリックフェーズの結果を用いたシミュレート画
像，ノンパラメトリックフェーズの結果を用いたシミュレート画像の画素 (u, v)
の輝度値をそれぞれ iobs (u, v), isim,p (u, v), isim,np (u, v) で表した場合に式 (45) で
示すように計算される誤差の改善値 iimp (u, v) を色相で示した．色相に関して，
緑の画素は誤差が変わらない部分で，赤色に近づくほど誤差の改善が顕著に見
られる画素となっている．この画像上での赤色や黄色の部分において，誤差評
価関数の値が減少したと考えられる．
49
図 58: 探索 ID1 の pars でのシミュレート画像
iimp = |iobs (u, v) − isim,p (u, v)| − |iobs (u, v) − isim,np (u, v)|
50
(45)
図 59: 探索 ID1 の pare でのシミュレート画像
51
図 60: 観測画像，pars でのシミュレート画像，pare でのシミュレート
画像の比較
図 61: ノンパラメトリックフェーズでの誤差の改善
52
第6章
おわりに
本研究ではチンダル現象下で実物体にレーザポインタを照射した様子を観測
することで，これまで反射光がカメラで観測できないために計測が困難であっ
た金属物体や黒色物体についても計測が可能となり，さらに，画像中で観測され
る入射光および反射光の振る舞いから物体表面の数式モデルを用いない BRDF
も同時に獲得する手法を提案した．
三次元形状に関しては提案手法により精度よく推定することができた．一方，
各反射方向への反射率及び散乱パラメータから生成したシミュレート画像と観
測画像の誤差評価関数を最小にするパラメータの探索を試みたが，誤差評価関
数がノイズに対して頑健でない設計であり，複数のパラメータを適度に変化さ
せるとそれらの変化の相殺によりほぼ同じシミュレート画像が生成されるとい
う原因で局所解に陥ってしまった．
今後の課題としては，誤差評価関数をノイズに頑健となるように設計するこ
とや各反射方向への反射率と散乱パラメータを別々に推定することなどにより，
各々のパラメータが一意に収束するようにすることが挙げられる．
53
謝辞
本研究を進めるにあたり，熱心なご指導を賜りました美濃導彦教授に深く感
謝いたします．本論文の副査として多大なアドバイスを賜りました中村裕一教
授に深く感謝いたします．また，研究の方向性について日頃よりアドバイスを
賜りました椋木雅之准教授，飯山将晃准教授，舩冨卓哉助教に深く感謝いたし
ます．最後に，グループミーティングでアドバイスをして下さいましたモデル
グループの皆様，ならびに美濃研究室の皆様に深く感謝いたします．
54
参考文献
[1] 大田 雄也, 川崎俊央, 佐川 立昌, 古川 亮, 川崎 洋, 八木 康史. デブルーイン
系列と Belief-Propagation を用いた高密度ラインパターン検出による高速
動体の 3 次元計測手法. 第 12 回画像の理解・認識シンポジウム (MIRU2009)
論文集, pp.222-229, Jul 2009.
[2] M. Lindstrom, J. -O. Eklundh. Detecting and tracking moving objects from
a mobile platform using a laser range scanner. Proc. of IEEE Int. Conf.
on Intelligent Robots and Systems, Vol.3, pp.1364-1369, 2001.
[3] 舩冨 卓哉, 飯山 将晃, 角所 考, 美濃 導彦. 光切断法による複数運動剛体
の形状計測. 電子情報通信学会論文誌, Vol.J90-D No.8 pp.1858-1867, Aug
2007.
[4] 町田 隆, 竹村 治雄, 横矢 直和. インバースレンダリングによる物体表面反
射特性の推定. 情報処理学会論文誌, Vol.47, No.SIG 9, pp.69-86, Jun 2006.
[5] Hongsong Li, Sing-Choong Foo, Kenneth E. Torrance, Stephen H. Westin.
Automated three-axis gonioreflectometer for computer graphics applications. Proc. SPIE, Vol.5878, pp.221-231, 2005.
[6] Yasuhiro Mukaigawa, Kohei Sumino, Yasushi Yagi. Rapid BRDF Measurement Using an Ellipsoidal Mirror and a Projector. IPSJ Transactions on
Computer Vision and Applications. Vol 1, pp.21-32, Jan 2009.
[7] 小俣 和子, 斎藤 英雄, 小沢 慎治. 光源の相対的回転による物体形状と表面
反射特性の推定. 電子情報通信学会論文誌, Vol.J83-D2, No.3, pp.927-937,
Mar 2000.
[8] K. E. Torrance, E. M. Sparrow. Theory of off-specular reflection from
roughened surfaces. Journal of the Optical Society of America. 57, pp.11041114, 1967.
[9] B.T. Phong．llumination for computer generated picture．Communications of the ACM, vol.18, no.6, pp.311-317, June 1975．
[10] Zhengyou Zhang．A Flexible New Technique for Camera Calibration．
IEEE Trans. PAMI, Vol.22, No.11, pp.1330-1334, 2000.
[11] http://opencv.jp/
[12] T. Ueshiba and F. Tomita．A Factorization Method for Projective and
55
Euclidean Reconstruction from Multiple Perspective Views via Iterative
Depth Estimation．5th European Conference on Computer Vision, pp.296310, 1998.
[13] Srinivasa G. Narasimhan, Mohit Gupta, Craig Donner, Ravi Ramamoorthi, Shree K. Nayar, Henrik Wann Jensen. Acquiring scattering properties of participating media by dilution. ACM SIGGRAPH. Vol.25, No.3,
pp.1003-1012, 2006.
[14] Norihisa Segawa, Yukio Yasuhara, Naohisa Sakamoto, Tomoki Yoshihisa,
Yasuo Ebara, Koji Koyamada. A Real-time Sensor Network Visualization
System using KVS - Kyoto Visualization System. the 5th ACM Conference
on Embedded Networked Sensor Systems (Sensys’07), pp.367-368, 2007.
[15] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery. Numerical Recipes the Art of Scientific Computing Third Edition.
[16] 馬場葉子, 向川康博, 八木康史. 照明光と反射光の散乱を考慮した反射特性
の推定. 信学技報 Vol.110 No.381 PP.237-244, 2011.
[17] http://www.betterseishin.co.jp/download/pdf/kussetu.pdf
56